Luis Cayuela

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

2012-02-01T23:09:00.000+01:00

No se ha hecho esperar la nueva versión1.1. del paquete TPL (anteriormente v. 1.0) para realizar laestandarización taxonómica de nombres de plantas utilizando unaconexión en línea al portal de The Plant List. Gracias a loscomentarios de varios colegas, he corregido unos cuantos bugs, afinado losmecanismos 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 últimono es de interés para el usuario estándar, pero para aquellos quequieran utilizar el código fuente para modificarlo conposterioridad, sí será de gran utilidad, ya que el código ahora esalgo más legible y está mejor estructurado.

He probado el paquete con varias basesde datos de plantas, una propia que he sacado del proyectoBIOTREE-NET, con más de 5000 nombres de árboles para todaCentroamé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 nombresde briófitos que me ha pasado Íñigo de la Cerda-Granzow, unalista de 238 árboles de la Reserva del Triunfo, en México, que me hapasado 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 paratodas ellas.

¡Ojo! Esto no quiere decir necesariamente que el códigono 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ódigofuente para ver la gran variedad de situaciones que podemosencontrarnos). Aunque he hecho comprobaciones manuales (remitiéndomea TPL) con nombres elegidos al azar y aparentemente todo funcionabien, sería conveniente que los usuarios finales haganverificaciones de este tipo, sobre todo en estas primeras etapas, yme comuniquen si encuentran alguna incoherencia entre los resulutadosobtenidos en R y lo que realmente dice TPL.

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

Otra de las mejoras de esta nuevaversión del paquete TPL es que el código corrige ahora con másexactitud los errores tipográficos en los nombres de las plantas.Los errores tipográficos son más comunes de lo que creemos. Esfácil poner una letra de más o de menos en un nombre complejo. Aveces, 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íaerrores tipográficos, y la función TPL permite corregir conbastante exactitud muchos de ellos (siempre y cuando los erroresestén en el epíteto específico). En el listado de más de 5000plantas del proyecto BIOTREE-NET, identificamos, por ejemplo, 299errores, que fueron corregidos automáticamente. Se pueden ver todosestos nombres y el consiguiente output de la función TPL aquí.Aparentemente todo está bien.

Bueno, creo que ya me he enrollado bastante. Elpaquete 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, cualquierincidencia, sugerencia o comentario hacédmelo llegar. Y si todofunciona bien, también agradecería comentarios al respecto (einformación sobre el número de nombres con el que se ha probado la función,etc).

¡Cruzad los dedos y... a correr elcódigo!

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

2012-01-27T13:12:00.000+01:00

Cuando se trabaja con grandes bases dedatos de vegetación con procedencia muy diversa, la taxonomía puedejugarnos una mala pasada. En estas bases de datos es frecuenteencontrar: (1) especies que se llaman de distinta forma pero que enrealidad son la misma (sinónimos); (2) especies con el mismo nombrepero que en realidad son distintas (homónimos); (3) errorestipográ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 entodo este caos que suponen los sistemas nomenclaturales (aunque obviamente no está exenta de errores).

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

Pues bien, en el marco del proyectoBIOTREE-NET, una iniciativa que trata de compilar y estandarizarinformación de árboles en inventarios forestales para todaCentroamérica, he diseñado una función en R que te permite hacertodo este trabajo de forma automatizada. Esto ahorra mucho trabajomanual y puedes cotejar miles de nombres en poco tiempo (unas horascomo mucho). La función (TPL) se presenta dentro de paquete de Rcon el mismo nombre, y utiliza a su vez otra función (TPLck) quehace el cotejo para nombres individuales.

El procedimiento está resumido en elsiguiente 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 plantases un arreglo de datos (data.frame) con información sobre el estatusdel nombre según The Plant List, el nombre actualizado (en caso desinónimos o errores tipográficos), la familia y la autoridad, entre otras cosas.

Las principales limitaciones hasta elmomento son que: (1) no es posible resolver el problema de loshomónimos; (2) las correcciones de errores tipográficos sólo serealizan cuando los errores existen en el epíteto específico. Silos errores se producen en el género o en el epítetoinfraespecífico, no hay nada que hacer; (3) en el caso de que unaespecie sea sinónima de otra y mantenga el mismo nombre cambiandosólo la autoridad (ej. Bartramia pomiformis var. elongata Turner comosinónimo de Bartramia pomiformis var. elongata Hedw.), el procedimientoextráe finalmente la autoridad del sinónimo (Turner) y no delnombre aceptado (Hedw.); (4) en el caso de que haya muchos epítetosinfraespecíficos y ninguno se corresponda con el nombre que estamosvalidando, la función busca por defecto el nombre de la especie queNO tenga epíteto infraespecífico. Si este nombre no existe en TPLel nombre se queda sin resolver (pero se da una advertencia alusuario para que pueda revisar este nombre a posteriori).

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

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

2012-01-16T13:34:00.000+01:00

El editorial de la revista científicaNature 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 beneficiosdesproporcionados. 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 Europamás competitiva.

Dicho y hecho; continuando la tradición de los últimos años delgobierno Zapatero, el nuevo Gobierno de Rajoy en España ha anunciadoun recorte de 600 millones de Euros en I+D+i. El presupuesto para2011 (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 recortedeja a la Ciencia española en una situación de emergencia, en uncontexto de crisis económica y política nacional. También hay quetener en cuenta al estudiar esas cifras que buena parte del I+Despañol financia investigación en tecnología militar, con lo que el dinero dedicado al tipo deinvestigación que querrían mayoritariamente los españoles esincluso menor.La investigación y la innovación son pilares fundamentales parael desarrollo de una sociedad moderna, y hemos demostrado muchasveces que los científicos españoles podemos ser tan buenos como losdel resto de países si nos dan la oportunidad. Por ello, elcientífico español Francisco J. Hernández ha impulsado unainiciativa 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)

2011-12-15T15:40:00.000+01:00

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

2011-11-29T22:48:00.000+01:00

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

2011-11-24T14:21:00.004+01:00

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?

2011-11-04T22:11:00.007+01:00

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

2011-10-26T09:45:00.006+02:00

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:

Agotan la fertilidad del suelo.
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.
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.
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

2011-10-13T19:14:00.007+02:00

¿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

2011-09-29T22:37:00.003+02:00

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

2011-08-17T11:39:00.008+02:00

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 festín de orugas

2011-07-06T12:51:00.008+02:00

La procesionaria del pino (Thaumetopoea pityocampa) es un lepidóptero típico de la región Mediterránea. La mariposa de la procesionaria se aparea en verano. La hembra pone sus huevos sobre las copas de los árboles y 30 ó 40 días después nacen las orugas, generalmente en los meses de agosto y septiembre, que construyen sus nidos sobre las ramas y pasan el invierno en ellos. Entre febrero y abril descienden al suelo, forman las características filas indias –de ahí su nombre común de “procesionaria”– y se entierran finalmente en el suelo, donde pasan a la fase de crisálida. En verano las crisálidas eclosionan y surgen las mariposas, que se aparean y reinician de nuevo el ciclo. Durante el invierno, las orugas se alimentan de las hojas de los pinos en los que construyen sus nidos y esa es precisamente la causa de la defoliación.

A pesar de su toxicidad, existen varias especies de aves que han desarrollado mecanismos y estrategias para poder alimentarse de las larvas de procesionaria. Así, por ejemplo, el críalo europeo (Clamator glandorius) y el cuco (Cuculus canorus) son capaces de regurgitar los pelos urticantes de las orugas. El herrerillo capuchino (Lophophanes cristatus) y el carbonero común (Parus major) y garrapino (P. ater) no son capaces de ingerir la oruga entera, sino que las pelan como pipas, quitándoles la cabeza y el tegumento para alimentarse solamente de la parte carnosa de la larva. Estos últimos son los principales responsables de los agujeros que vemos en los nidos de la procesionaria. Existen otras especies que se alimentan de la procesionaria, pero no durante la fase de oruga, como la abubilla (Upupa epops; ver foto más abajo) que desentierra las crisálidas del suelo, o el chotacabras gris (Caprimulgus europaeus), que es capaz de cazar a la mariposa al vuelo durante su efímera existencia (generalmente no más de 24 horas).

Todas estas especies actúan como agentes efectivos para el control natural de este insecto, evitando así con su presencia la aparición de superpoblaciones de procesionaria en los pinares mediterráneos.

Barbaro, L., & Battisti, A. (2011). Birds as predators of the pine processionary moth (Lepidoptera: Notodontidae) Biological Control, 56, 107-114 : doi:10.1016/j.biocontrol.2010.10.009

La Fundación Gypaetus: Un ejemplo de mala gestión

2011-06-30T15:14:00.006+02:00

La Fundación Gypaetus que desde hace años realiza de forma exitosa programas de reintroducción del quebrantahuesos, sensibilización social y seguimiento y lucha contra el veneno en Andalucía, se ve ahora amenazada por la mala administración de su gerente, Jesús Charco, quien ha llevado a cabo una serie de despidos en masa no justificados, frenando muchos de los programas de conservación que la Fundación venía desarrollando. Las decisiones tomadas no solo han afectado a los trabajadores de la Fundación, sino que han repercutido negativamente en fondos públicos, y lo peor de todo, han acabado afectando al éxito del proyecto de reintroducción del quebrantahuesos en Andalucía, de manera que en la prensa leemos que “este año no habrá suelta de quebrantahuesos para evitar envenenamientos”.

Despidos en masa

En solo tres meses la Fundación Gypaetus ha prescindido de todo el equipo técnico de la sede central de Jaén: 6 personas han sido despedidas y a 3 no les han renovado el contrato. Estos despidos en bloque (todos con reconocimiento de improcedencia, lo que conlleva la respectiva indemnización) afectan a personal técnico que llevaba años trabajando en programas de liberación del quebrantahuesos, sensibilización, contacto y fortalecimiento de relaciones con los agentes sociales involucrados en las áreas, como ganaderos y cazadores. Esos despidos suceden inmediatamente después de que este equipo hubiese escrito una carta al Patronato de la Fundación Gypaetus mostrando su preocupación ante la gestión del nuevo gerente, Jesús Charco, que había derivado en insostenibles problemas económicos y de otra índole a partir de su ingreso en la Fundación a finales de 2009. Ese fue el mismo equipo que logró que la Comisión Europea financie durante cinco años (iniciándose en 2010) el proyecto que idearon, presentaron y consiguieron, “Innovative actions against ilegal poisoning in EU Mediterranean pilot areas”, y el mismo equipo que no hace mucho recibió por parte de la Comisión Europea el reconocimiento al proyecto Life ejecutado entre 2004 y 2009, considerándolo entre los "Best Life Projects 2010".

No es comprensible que una Fundación considere necesario deshacerse de personas que, no sólo ganan premios a nivel europeo, sino que llevan a la institución millones de euros para que puedan seguir ejecutándose sus programas.

La indemnización que se ofrecía tras los despidos, en cada uno de los casos, fue inferior a la que la legislación reconoce, lo que ha originado procesos judiciales, que han desembocado en el reconocimiento por sentencia de la indemnización reclamada y en el pago por parte de la Fundación Gypaetus de salarios de tramitación. Esto supone que a lo largo de los últimos meses dicha Fundación ha desembolsado casi 90.000 € para cubrir estos despidos improcedentes, lo que supone un 16% del presupuesto anual para la reintroducción del quebrantahuesos. Estos despidos se han pagado con dinero público destinado a programas de reintroducción del quebrantahuesos, lo que ha repercutido negativamente en el desarrollo de dichos programas. Además de paralizar las liberaciones de pollos y de suprimir las tareas de hacking, las cuatro personas que trabajaban en el programa de seguimiento y lucha contra el veneno han sido sustituidas por una sola, a la que además se le han dado instrucciones de no salir al campo ¿cómo vigilar y evitar así las posibles amenazas? Esto ha contribuido al envenenamiento de dos quebrantahuesos en la Sierra del Castril. La denuncia de esta situación le ha costado su puesto a uno de los nuevos trabajadores en mayo de 2011.

Como simular un bosque en 3D

2011-03-15T14:29:00.004+01:00

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

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

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

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

Created by Pretty R at inside-R.org

Y el resultado sería este:

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

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

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

2011-01-25T23:35:00.011+01:00

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

El NMDS se implementa de la siguiente forma:

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

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

Problemas de convergencia

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

La persistencia de la lactasa en Europa

2010-12-14T15:45:00.005+01:00

La mayoría de las personas adultas en el mundo no producen la enzima lactasa y por tanto no son capaces de digerir la lactosa que se encuentra en la leche. Esta enzima está presente en los recién nacidos, pero desaparece a la edad del destete, entre los dos y los tres años, edad en la que ya no se hace necesario digerir la lactosa que se encuentra presente en la leche materna. En Europa, así como en algunas poblaciones pastorales de África, Oriente Medio y el sur de Asia, se da la condición de que muchos adultos siguen produciendo lactasa. A este fenómeno se le conoce como "persistencia de la lactasa". En Europa esta condición se produce entre el 15% y el 54% de la población de los países del sur y este de Europa, entre el 62% y el 86% de la población de Europa central y occidental y entre el 89% y el 96% de la población de las Islas Británicas y Escandinavia. La persistencia de la lactasa está asociada con una mutación en el gen MCM6 a la que se la denomina -13,910*T. Curiosamente, la persistencia de la lactasa que manifiestan las poblaciones pastorales de África, Oriente Medio y el sur Asia no está asociada a este gen, sino a otros, lo que apunta claramente a un fenómeno de convergencia evolutiva.

Este estudio investiga mediante técnicas genéticas, datos arqueológicos y sofisticadas técnicas de análisis estadístico el origen de la persistencia de la lactasa en Europa. Los resultados muestran que la aparición de la mutación -13,910*T apareció hace aproximadamente unos 7500 años y se extendió como consecuencia de procesos de selección natural. Dado que el consumo de leche fresca no fue posible sin la domesticación de los animales, es muy probable que la persistencia de la lactasa coevolucionara con la práctica cultural del ordeño. Aunque el epicentro geográfico de la mutación se encuentra actualmente en el norte de Europa, los resultados del estudio parecen indicar que se originó en la región de los Balcanes y Europa central y que posteriormente se extendió hacia el norte de Europa.

Referencia

Itan Y, Powell A, Beaumont MA, Burger J, & Thomas MG (2009). The origins of lactase persistence in Europe. PLoS computational biology, 5 (8) PMID: 19714206

¿Afecta lo que comes al medio ambiente?

2010-10-19T11:09:00.007+02:00

Más allá del tópico "¿eres lo que comes?" nuestra dieta puede influir no sólo a nuestra salud, sino también al medio ambiente. Varios estudios científicos dan evidencias de ello.

El estudio de Luigi Baroni (Baroni et al. 2007) realizado con base en distintos tipos de dietas y sistemas de producción en Italia demuestra que las dietas omnívoras tienen un impacto mucho mayor que las dietas vegetarianas o veganas. La razón de ello radica en el hecho de que es mucho más costoso (energética y ambientalmente) producir una caloría de carne (que requiere 40 calorías de combustible) que una caloría de vegetal (que requiere en promedio 2.2 calorías de combustible). La producción de carne (y en menor medida también la producción de otros productos derivados como la leche o el queso) es, por tanto, muy ineficiente y consume una gran cantidad de recursos naturales, entre ellos combustibles fósiles. En Europa, por ejemplo, producimos proteínas vegetales suficientes para alimentar a toda la población europea, pero no tenemos proteínas vegetales suficientes para alimentar a todo el ganado, por lo que cerca del 80% de los productos vegetales utilizados en alimentación animal se importan de otros países. En estos países la producción intensiva de cultivos (como la soja en Brasil y Argentina) ha reemplazado a muchos sistemas extensivos de policultivos o, lo que es mucho más alarmante, bosques nativos como la selva Amazónica, cuya destrucción es causada en gran medida por la producción de grano y forraje para la alimentación del ganado de Estados Unidos y Europa.

Uno de los mayores impactos en lo que a consumo de recursos naturales se refiere es el consumo de agua. Las actividades agropecuarias son responsables del 70% del consumo de agua dulce en todo el planeta, mientras que sólo el 22% se utiliza en la industria y el 8% se destina a uso doméstico. La mayor parte del agua utilizada en estas actividades se destina a la producción (a veces intensificada mediante la implementación de sistemas de regadío) de cereales o plantas oleaginosas (soja, girasol, algodón, lino, etc.), de las cuales cerca del 50% son destinadas al consumo de ganado. Una gran cantidad de agua es también consumida directamente por el ganado y, una parte nada despreciable se dedica al mantenimiento (limpieza) de establos, centrales lecheras, mataderos, etc. ¿Y qué hacemos con los excedentes que, en forma de excretas, se derivan de los sistemas intensivos de producción de carne? En un sistema de producción extensivo, estas excretas suponen un fertilizante natural muy potente, pero cuando el ganado está estabulado, los residuos generalmente crean un problema ambiental importante, pudiendo llegar a contaminar el agua de los ríos e incluso de los acuíferos.

Otro estudio más actual (Deckers 2010) apunta a que el principal impacto derivado de los sistemas actuales de producción de carne sería la emisión de gases de efecto invernadero a la atmósfera, lo que contribuiría en cerca de un 50% al calentamiento global. De seguir creciendo la población mundial y el consumo de carne, para el 2050 la contribución de la ganadería a las emisiones de gases de efecto invernadero sería de entre el 66.9% y el 83.7% según distintos escenarios. Y esto asumiendo de forma optimista que la producción de carne generaría un 20% menos de emisiones como consecuencia del aumento en eficiencia de los sistemas de producción. En otro artículo publicado en la prestigiosa revista The Lancet, Anthony J. McMichael y colaboradores (2007) proponen que para reducir las emisiones de gases de efecto invernadero de aquí al año 2050 se debería reducir el consumo promedio de carne per cápita de 100 gr a 90 gr diarios. Parece una reducción muy pequeña, pero hay que considerar que en la mayoría de los países industrializados, el consumo promedio de carne por día es muy superior a esta cifra (en el caso de Estados Unidos por ejemplo, el consumo promedio por día es de 218 gr).

Otros argumentos a favor de disminuir el consumo de carne tienen que ver con la salud, incluyendo aquí el uso de sustancias nocivas como antibióticos y hormonas que se le suministra al ganado estabulado en sistemas industrializados de producción (Sapkota et al. 2007). Algunas personas también utilizan argumentos éticos en relación al sufrimiento de animales bajo los sistemas de producción intensivos. A estos y otros respectos, el libro recientemente publicado de Joyce D'Silva y John Webster (2010) se convierte una referencia imprescindible.

Resumiendo: el principal problema no es consumir carne, sino que los sistemas de producción de carne son más inefectivos que los de producción vegetal. Si a esto sumamos que la población mundial sigue creciendo y que el consumo de carne per cápita también crece (se ha multiplicado por cuatro en el último medio siglo), el problema se hace más que patente. Nos movemos (en realidad ya estamos) hacia un modelo insostenible de producción, por lo que hace falta promover un cambio y generar nuevas alternativas. ¿El cambio? Disminuir nuestro consumo de carne a nivel individual ¿Las alternativas? Sistemas de producción orgánicos que tengan un menor impacto ambiental.

Referencias

Revistas relacionadas con Producción de Materias Primas

2010-09-08T09:32:00.003+02:00

Este año he empezado a impartir una asignatura nueva en el nuevo grado de Ciencia y Tecnología de los Alimentos titulada Producción de Materias Primas. Como este es un tema relativamente nuevo para mi, estoy investigando qué revistas científicas publican cosas interesantes en el área de agricultura y alimentación. He aquí un primer listado (se agradecerán nuevas aportaciones) de revistas que pueden ser de interés en la materia:

También es muy interesante el Food Reviews, publicado por el Economic Research Service del Departamento de Agricultura de los Estados Unidos.

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

2010-09-07T14:39:00.011+02:00

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:

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);
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);
Se ordenan estas áreas en orden descendente, siendo las más importantes la que más riqueza, rareza o vulnerabilidad tienen.
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:

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

¿Es la leche sana y saludable?

2010-08-03T14:26:00.023+02:00

En la literatura científica existe, desde hace años, un debate intenso sobre el papel de la leche como fuente principal de calcio en la dieta occidental. Muchos estudios cuestionan este modelo y algunos, incluso, van más allá, demostrando que el consumo de productos lácteos es desaconsejable debido a la intolerancia a la lactosa, común en muchas personas, y puede favorecer la aparición de diversas enfermedades, como la diabetes A. Este debate, por desgracia, no ha permeado al resto de la sociedad, que sigue creyendo en un modelo único que predica que la leche es buena y que, sin la leche, no tenemos suficiente cantidad de calcio en nuestra dieta, sin cuestionar ni dudar por un solo momento de su validez.

Cuando por primera vez escuché que la leche no era buena para la salud mi primera reacción fue la negación. Mi negación se amparaba en el inconsciente colectivo, como muchas otras supuestas verdades o verdades cuestionables, que forman parte de nuestra manera de ser, como individuos y como sociedad. Yo, como la mayoría de la gente de mi generación (y de generaciones anteriores y posteriores), me he criado tomando leche de vaca, y por tanto cuestionar esta verdad supone asumir un error sistemático en nuestra forma de alimentación. Asumir que nos hemos equivocado -o que podemos estar equivocándonos- es mucho más difícil todavía cuando se trata de la salud de nuestros hijos. Por eso mucha gente que da leche de vaca a sus hijos no quiere siquiera hacerse esta pregunta.

Cuando por segunda vez escuché esta misma idea, mi reacción fue negarla otra vez. Mi principal argumento para desmontar a mi interlocutor fue pedirle que me dijera si estas afirmaciones estaban basadas en trabajos científicos. Claro, con esto tenía la batalla ganada, porque la mayoría de la gente carece de esta información. En mi interior, no obstante, se había sembrado la semilla de la duda y empecé a cuestionar si realmente la creencia bien establecida de que la leche es un alimento bueno y casi necesario para un correcta alimentación tenía fundamento. Cuando nació mi hijo hace algo más de un año, leí algunos libros de pediatría (ej. Gutman 2009) que me acabaron de convencer de los efectos perniciosos de la leche. Los argumentos que se daban en estos libros en contra del consumo de productos lácteos estaban supuestamente respaldados por estudios científicos (aunque no se especificaba de manera explícita). Consecuentemente mi mujer y yo decidimos que nuestro hijo no tomaría leche de vaca de manera habitual.

Como científico (aunque de otra disciplina, la ecología) soy consciente de las limitaciones de la aplicación del método científico para extraer conclusiones unívocas en relación al objeto de estudio. A veces dos estudios científicos muestran resultados totalmente contrarios con respecto a una misma pregunta. Otras, los resultados son interpretados libremente más allá del ámbito del estudio para apoyar determinadas posturas o acciones por parte de políticos, activistas, empresarios, etc. Por ello hoy, casi cuatro años después de escuchar por vez primera que el consumo de la leche no es tan bueno como nos hacen creer, he decidido buscar las evidencias científicas que respaldan una u otra teoría.

Lo más sorprendente de mi búsqueda es que la información es muy accesible, y por eso todavía me pregunto como el debate no se ha movido al ámbito social y por qué mucha gente todavía se sorprende cuando oye esto por primera vez. Me ceñiré aquí al ámbito científico de mi búsqueda. En cuanto entras en Google académico y escribes algunas palabras clave como "Cow milk" + "diet" aparecen un gran número de referencias a artículos científicos publicados en prestigiosas revistas médicas como Pediatrics, American Journal of Public Health o The American Journal of Clinical Nutrition, que hablan de la intolerancia a la lactosa, la absorción de calcio en la leche frente al calcio absorbido en otros alimentos (de origen vegetal), la leche como causante del constipado crónico, cáncer de próstata, etc. A no ser que se tenga una suscripción a estas revistas, el acceso a estos artículos no es posible, pero sí se puede consultar el resumen del artículo, que contiene una descripción de los métodos usados y de las principales conclusiones. Aún así, el lenguaje es difícilmente accesible para quién no esté acostumbrado a la jerga científica. Mi intención aquí es hacer un resumen de los principales aspectos relacionados con el consumo de la leche aportando evidencias científicas y con referencias específicas a (algunos de) los estudios en los que se basan estas afirmaciones. Esta información, como he dicho antes, no es nueva. Se pueden encontrar revisiones del tema en castellano e inglés en numerosas páginas de internet, entre las que destacaré la de José Ramón Llorente (Presidente de la Asociación Española de Nutrición Ortomolecular) "La leche, ese producto pernicioso para los seres humanos" y la página del School of Public Health de la Universidad de Harvard "Calcium and milk: what's best for your bones and health?".

La leche y el calcio

La razón fundamental por la que los nutricionistas occidentales (no así los orientales) recomiendan tomar leche y sus derivados es porque la consideran muy nutritiva y especialmente rica en calcio, agregando que la ingesta periódica de ese mineral es imprescindible para mantener la salud, sobre todo la de los huesos. La confrontación de estudios a favor y en contra de la leche radica en la cantidad requerida de calcio y en su disponibilidad para ser absorbida por el organismo. Así, en EEUU la cantidad diaria recomendada de calcio para niños es de > 1300 mg. En España, el RD 1669/2009 indica que la cantidad diaria recomendada (CDR) de calcio debe ser de 800 mg, sin hacer distinción entre niños y adultos. No todo el calcio que ingerimos con los alimentos está disponible para el organismo. Varios autores (Heaney et al. 1988, Charles 1992, Weaver y Pawecki 1994, Weaver et al. 1999) coinciden en que la leche y, en general los productos lácteos, no sólo tienen mucho calcio, sino que además este calcio está disponible en una gran proporción (alrededor de un 30%) en comparación con otros alimentos como las espinacas, el brócoli, algunos cereales, el tofu o la leche de soja, que son ricos en calcio, pero con una baja disponibilidad para el organismo (ver revisión en Goldberg et al. 2002). Muchos estudios relacionan el consumo de leche en edades tempranas con una disminución en el riesgo de osteoporosis en la edad adulta (ver una revisión en Heaney 2000). Otros estudios, por el contrario, apuntan a que ni la leche ni la ingesta de calcio por encima de los 500 mg diarios parecen tener un efecto beneficioso sobre los huesos de los niños (ej. Feskanich et al. 1997, Lloyd et al. 2000, Kanis et al. 2005). Una revisión de 37 estudios (Lanou et al. 2005) en dónde se relacionaba el consumo de lácteos con la salud de los huesos evidenció una ausencia de relación entre la ingesta de leche y varios indicadores de la salud de los huesos en la mayoría de los casos (28), mientras que, en unos pocos casos (9), se detectó un efecto positivo (si bien muy pequeño) del consumo de leche en la salud de los huesos.

Lo que debemos de comprender es que los riesgos de fractura de huesos no están necesariamente ligados con la ingesta de calcio. Unas cantidades mínimas de calcio son requeridas, por supuesto, pero hay otros factores que pueden ser más importantes en la absorción del mismo. Uno es el ejercicio físico que, en dosis moderadas, favorece notablemente la absorción de calcio (Lloyd et al. 2000). Otro parece ser la ingesta de proteínas. Algunos estudios han relacionado las dietas altas en proteína de origen animal y el sodio (sal) con una descalcificación de los huesos (Barzel y Massey 1998, Frasseto et al. 2001). La Organización Mundial de la Salud (OMS), explica que el consumo elevado de proteínas, sobre todo de origen animal, podría contrarrestar los efectos de una alta ingesta de calcio en la dieta" y recomienda el ejercicio físico, reducir la ingesta de sales y aumentar el consumo de frutas y verduras para favorecer el fortalecimiento de los huesos (OMS/FAO 1998). Así se da la paradoja de que la mayor incidencia de osteoporosis se produce en EEUU y en países occidentales, que son los principales consumidores de leche a nivel mundial, mientras que permanece bajo en Asia y África donde el consumo de leche es realmente bajo o, en muchos casos, inexistente.

En conclusión, no está muy claro que la ingesta de leche realmente satisfaga las necesidades de calcio debido al problema de su absorción. En países con dietas muy ricas en proteínas animales y sales, como EEUU, los requerimientos de calcio serán mayores que en países con dietas más equilibradas y ricas en frutas y verduras. Si la leche no fuera asociada a determinadas enfermedades, como veremos más adelante, no habría mayor problema en cuanto a obtener la fuente de calcio de la leche o de otros alimentos, lo cual está comprobado que es posible. Sin embargo, la polémica no acaba aquí.

La intolerancia a la lactosa

Los bebés tienen una enzima, que es la lactasa, que permite metabolizar el azúcar de la leche (la lactosa). Este enzima va desapareciendo con el destete, más o menos a partir del año, aunque parece que en la raza blanca permanece por más tiempo que en la raza negra o asiática (Scrimshaw y Murray 1988). Cuando carecemos de estas enzimas, la lactosa deja de ser metabolizada en el intestino delgado y pasa al intestino grueso, dónde la lactasa es atacada por bacterias, produciendo fermentaciones (gases), cólicos, diarreas, etc.

La intolerancia a la lactosa es la norma más que la excepción (Barnard 2003). Muchos problemas de digestión o cólicos frecuentes podrían tener su origen en el consumo de lácteos. En otros casos, la respuesta no es tan inmediata ni tan visible, y las personas que tienen intolerancia a la lactosa pueden tolerar relativamente bien los productos lácteos (Goldberg et al. 2002), sin que se produzcan ninguno de estos síntomas a corto o medio plazo. Sin embargo, cabe la duda de si pueden aparecer efectos a largo plazo. Por otro lado, algunas personas sí tienen lactasa y son por tanto capaces de digerir este carbohidrato, por lo que para ellos el consumo de leche no sería ningún problema. Esto es más frecuente en las poblaciones del norte de Europa, mientras que en los países mediterráneos se estima que cerca del 90% de la población tiene intolerancia a la lactosa.

Mucosidades y alergias

Además de la lactosa, que es algo que la leche de vaca comparte con la leche humana, la leche de vaca tiene un 300% más de caseína que la leche humana. La caseína es una proteína muy densa que obstruye el sistema respiratorio. La reacción más común de nuestro sistema inmunológico frente a la absorción de proteínas extrañas es la secreción de moco en la nariz y la faringe. La acumulación constante de mocos puede agravar un resfriado común, derivando en rinitis, sinusitis, bronquitis, otitis, neumonía e infecciones de oídos. Diversos estudios han demostrado la relación entre el consumo de leche y el constipado crónico (Iacono et al. 1995) el reflujo gastro-intestinal (Salvatore y Vandenplas 2010), diversos tipos de alergia y otitis (Juntti et al. 1999, Rona et al. 2007).

Otras consideraciones

El segundo informe de expertos del World Cancer Research Fund y el American Institute for Cancer Research (WCRF/AICR 2007) sugiere que hay un mayor riesgo de sufrir cáncer de próstata en dietas con un alto contenido en calcio (como la leche, aunque no se refieren específicamente a ésta). También se ha establecido una relación entre el consumo de lácteos y la diabetes de tipo A (ej. Scott 1995, Virtanen et al. 2000), aunque algunos autores como Goldberg et al. (2002) advierten de posibles sesgos en estos estudios. Además de estos problemas, muchos vegetarianos y omnívoros están limitando su consumo de leche como consecuencia del aumento de sustancias tóxicas en la leche, como hormonas, antibióticos y pesticidas o, en un plano más ideológico, porque la producción de leche promueve una degradación ambiental (Steinfeld et al. 2006).

Algunos autores pueden llegar a ser muy alarmistas y relacionan el consumo de leche casi con cualquier enfermedad, como por ejemplo el autismo, pero la evidencia científica que respalda estos argumentos es, en la mayoría de los casos, poco sólida o inexistente.

Conclusiones

De toda esta polémica científica se pueden extraer algunas conclusiones:

La leche no es imprescindible para la ingesta diaria de calcio. En muchas sociedades no se consume leche y no tienen ningún problema en alcanzar las cantidades diarias recomendadas de calcio, alcanzando un nivel de incidencia de enfermedades relacionadas con los huesos, como la osteoporosis, mucho menor que en países dónde el consumo de lácteos es alto, como EEUU o muchos de los países europeos.
La importancia del calcio radica en su absorción, no en la cantidad ingerida, viéndose favorecida esta absorción por el ejercicio físico, las dietas bajas en sal y proteínas animales y ricas en frutas y verduras.
Una gran parte de la población tiene intolerancia a la lactosa, lo que puede manifestarse en dolores intestinales, gases, cólicos, etc., aunque hay gente con intolerancia a la lactosa que, aunque no es capaz de metabolizar la lactosa, no sufre ninguno de estos síntomas.
El alto contenido de la leche en proteínas como la caseína, en comparación con la leche humana, puede provocar mucosidad, catarros continuos, otitis y diferentes tipos de alergia.
La relación del consumo de leche con otras enfermedades como la diabetes A o el cáncer de próstata no está todavía muy clara, por lo que tenemos que ser cautos con como manejamos esta información sin caer en el alarmismo.

Como consumidores, no tenemos por qué dejar de lado el consumo de los productos lácteos. En lo que a mi respecta seguiré tomando queso y yogur, aunque posiblemente limite un poco más mi consumo a partir de ahora. El mayor problema lo veo en el grado de desinformación que tienen los supuestos expertos en salud (nutricionistas, dietistas, pediatras, etc.), que promueven alegremente el consumo de productos lácteos como parte de nuestra dieta básica. Esto es especialmente relevante en el caso de los niños, a quienes se recomienda el consumo de dos y tres vasos de leche al día, además de queso, yogurt, etc. ¿Es que estos profesionales no se actualizan, no leen, no tienen una mentalidad crítica ante los conocimientos aprendidos durante su época de formación, no se cuestionan que puedan estarse equivocando? Dar a nuestros hijos tres vasos de leche al día, solamente porque creemos que sin la leche no van a tener un crecimiento adecuado ni van a tener los huesos fuertes, es un acto sin fundamento, y en parte arriesgado, a la vista de las evidencias científicas. Debemos, eso sí, promover el consumo de otros alimentos ricos en calcio y promover una serie de pautas que favorezcan la absorción del mismo. Ahí lo dejo.

Análisis de datos ecológicos en R

2010-06-28T17:24:00.005+02:00

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

Los manuales pueden descargarse aquí:

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

Cursos de análisis de datos en ecología

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

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

BIOTREE-NET: Una red de inventarios forestales para la conservación de la biodivesidad en Centroamérica

2010-05-11T09:12:00.006+02:00

Ya se ha estrenado la página web del proyecto BIOTREE-NET, cuyo objetivo principal es consolidar una red de inventarios forestales en Centroamérica, con posibilidades de expansión hacia Sudamérica. La motivación de crear esta red es facilitar el intercambio de información entre científicos, gestores y conservacionistas con el fin de promover la investigación y la conservación de la biodiversidad en una de las regiones más diversas y amenazadas del planeta.

En el proyecto participan actualmente investigadores de diversas instituciones de Latinoamérica, Norteamérica y Europa, y se cuenta con datos de más de 1800 inventarios forestales distribuidos por el sur de México, Belize, Honduras, El Salvador, Nicaragua, Costa Rica y Panamá. Todos estos datos estarán disponibles en el futuro a través de un portal de datos que permitirá el acceso directo a la información biológica a través de una base de datos geográfica relacional.

Si estás interesado en el proyecto consulta www.biotreenet.com.

La magia de Gigapedia

2010-05-05T12:40:00.006+02:00

Hace poco un colega me enseñó como encontrar una gran cantidad de libros disponibles en formato electrónico. La respuesta estaba en Gigapedia: un portal de acceso a libros que cuenta con más de 40.000 libros disponibles. Gigapedia no mantiene los libros sino que da acceso a los mismos a través de enlaces ya existentes. Para poder hacer uso de esta valiosísima fuente de información y documentación es necesario registrarse en www.gigapedia.com. Una vez que nos hemos registrado, nos pedirán autentificar nuestro registro para lo cual tendremos que irnos a la dirección de correo que hemos suministrado durante el registro y activar la cuenta. Hecho esto ya tenemos acceso a Gigapedia y todos sus recursos.

Para hacer una búsqueda es importante que en el recuadro de búsqueda (situado en la parte superior derecha de la pantalla) seleccionemos Gigapedia. Sino por defecto hace una búsqueda normal en google.

Algunos libros interesantes

En Gigapedia he encontrado numerosos libros de ecología, biogeografía y bioestadística. Aquí van algunas referencias de libros de estadística (que utilizo para mis cursos) disponibles en Gigapedia:

Crawley, M.J. 2007. The R Book. John Wiley & Sons, UK.
Zuur, A.F., Ieno, E.N. & Smith, G.M. 2007. Analysing ecological data. Springer, New York.
Zuur, A.F., Ieno, E.N., Walker, N.J., Saveliev, A.A. & Smith, G.M. 2009. Mixed effects models and extensions in ecology with R. Springer, New York.

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

2010-04-14T11:45:00.007+02:00

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

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

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

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

Test de equivalencia ó cómo controlar el error de tipo II en inferencia estadística

2009-11-25T16:34:00.006+01:00

Muchas decisiones críticas para la conservación y gestión del medio natural se toman basándose en la creencia de que una diferencia estadística no significativa entre grupos significa que los grupos son iguales. Algunos experimentos se diseñan para ver si hay diferencias, por ejemplo, en el uso de dos estrategias de gestión para el control del fuego en un área protegida o para saber si dos poblaciones son genéticamente distintas entre sí. Los investigadores que llevan a cabo estos experimentos a menudo sacan conclusiones inapropiadas cuando no detectan diferencias estadísticamente significativas entre las poblaciones o grupos comparados, lo cuál puede tener repercusiones muy serias para la gestión y la conservación.

Los test de hipótesis clásicos normalmente se estructuran en torno a dos hipotesis: una “hipótesis nula” (representado por H0) que asume que no hay diferencias entre los grupos analizados, y una “hipótesis alternativa” (representada por H1) que generalmente establece que hay una diferencia “detectable” entre los grupos. Si la evidencia para H1 no es lo suficientemente fuerte, entonces se dice que H0 no puede ser rechazada. Sin embargo, muchos investigadores interpretan la falta de evidencia para rechazar H0 como evidencia para aceptar que los grupos comparados son iguales, lo cual es claramente incorrecto. ¿Por qué ocurre ésto? Cuando rechazamos la hipótesis nula, generalmente establecemos un valor de referencia α que nos indica el error de rechazar la hipótesis nula cuando no hay realmente diferencias entre los grupos (falso positivo o error de Tipo I). Por ejemplo, si α = 0.05 esto indica que 1 de cada 20 veces tendremos un falso positivo o, lo que es lo mismo, que podemos estar seguros en un 95% de que estamos en lo cierto cuando decimos que dos grupos son distintos entre sí. El problema viene cuando no rechazamos la hipótesis nula. En este caso, existen dos posibilidades: una es que realmente no haya diferencias significativas entre los grupos, y la otra es que haya diferencias entre los grupos pero que no la hayamos detectado (falso negativo), lo que se conoce como error de Tipo II o β. Por desgracia, mientras que para el error de Tipo I tenemos un cierto control de las probabilidades de equivocarnos, con los test de hipótesis clásicos, no hay manera de saber si al no rechazar la hipótesis nula estamos cometiendo un falso negativo o no.

Aunque el problema parece muy obvio, sigue habiendo un desconocimiento muy grande de la interpretación de los test estadísticos en el campo de la ecología y la biología de la conservación (y muy posiblemente en otros campos relacionados con la biología y las ciencias naturales). Una revisión de estudios publicados en las revistas Conservation Biology y Biological Conservation en 2003 encontró que casi dos tercios de las publicaciones con resultados no significativos interpretaron de manera inapropiada estos resultados como evidencia para decir que los grupos comparados eran homogéneos (Fidler et al. 2006).

¿Qué se puede hacer cuando un test estadístico no es significativo?

Algunos estudios sugieren el uso del poder estadístico (power analysis) para determinar el error de Tipo II (esto es, la probabilidad de equivocarnos al no rechazar la hipótesis nula). Sin embargo, muchos ecólogos no son conscientes de los problemas que estos test post-hoc plantean. Uno de estos problemas es que el poder observado está directamente (y negativamente) relacionado con el p-valor del test. Por lo tanto, cuando el p-valor no es significativo (p > 0.05) el poder observado será bajo. Por el contrario, cuando el p-valor sea significativo (p <0.05) style="font-weight: bold;">test de equivalencia -esto es, establecer como hipótesis nula que los dos grupos son diferentes y como alternativa que son iguales (Brosi & Biber 2009). El principal reto aquí es determinar qué diferencia mínima (Δ) entre los grupos es asumida como hipótesis nula. Si bien se puede pensar que este requerimiento introduce un elemento subjetivo en el análisis, también obliga al investigador a hacer sus supuestos sobre el proceso observado más explícitos. Por ejemplo ¿cuál es la mínima diferencia genética que un investigador puede asumir para determinar si dos poblaciones de una especie en peligro de extinción son distintas o iguales? Mediante un test de hipótesis tradicional podríamos llegar a la conclusión de que dos poblaciones son distintas, incluso aunque estas diferencias no tuvieran un significado relevante desde el punto de vista genético. Con el test de equivalencia formalizamos de alguna manera estas diferencias que nosotros, como investigadores, asumimos como relevantes desde el punto de vista biológico y no sólo estadístico, y además, conseguimos reducir el error de Tipo II (cómo se define en los test de hipótesis tradicionales) al 5% o incluso menos.

¿Cómo implementamos el test de equivalencia?

En realidad el test de equivalencia para la comparación de dos poblaciones no es más que un test de la t en dónde fijamos el parámetro mu (diferencia entre medias). En R la función tost() del paquete equivalence (Robinson 2008) permite hacer esto mismo, pero en el fondo podríamos llegar al mismo resultado utilizando la función t.test() del paquete stats. En la función tost() hay que especificar la diferencia mínima detectable, Δ, entre grupos, mientras que en la función t.test() tendríamos que definir el argumento mu y especificar como hipótesis alternativa una diferencia menor que la especificada (es decir, homogeneidad desde el punto de vista biológico). Esto último se especifica mediante el argumento alternative = “less”.

De igual modo suele ser bastante representativo ilustrar los intervalos de confianza de la diferencia en las medias. Esto nos da una idea de si los grupos son estadísticamente diferentes (intervalo no corta el cero), pero también de si los grupos son o no estadísticamente homogéneos de acuerdo a nuestro umbral mínimo detectable (intervalo queda comprendido entre ± Δ). Tomemos como ejemplo la figura de abajo (reproducida a partir de Brosi & Biber 2009) para ilustrar las posibles opciones. Podría ocurrir que:

(A, B) los grupos son diferentes (significación en el test de hipótesis tradicional) y además no homogéneos (no significación en el test de equivalencia);
(C) los grupos no son significativamente distintos (no significación en el test de hipótesis tradicional) pero son significativamente homogéneos (significación en el test de equivalencia);
(D) los grupos son significativamente distintos (significación en el test de hipótesis tradicional) y además son significativamente homogéneos (significación en el test de equivalencia. Esto puede ocurrir cuando las diferencias son detectables estadísticamente pero no relevantes desde el punto de vista biológico;
(E y F) los grupos no son significativamente distintos pero tampoco son significativamente homogéneos. Esto indica que son necesarios más datos para poder obtener una conclusión válida.

Berry, J. Brosi, & Eric G. Biber (2009). Statistical inference, Type II error, and decision making under the US Endangered Species Act Frontiers in Ecology and the Environment (7(9)), 487-494 : 10.1890/080003

Fidler, F., Burgman, M., Cumming, G., Buttrose, R., & Thomason, N. (2006). Impact of criticism of null-hypothesis significance testing on statistical reporting practices in conservation biology Conservation Biology, 20 (5), 1539-1544 DOI: 10.1111/j.1523-1739.2006.00525.x

A. Robinson (2008). Equivalence: provides tests and graphics for assessing test of equivalence. Package for the R Statistical Computing Language

Better species distribution modeling needed for the tropics

2009-10-05T19:39:00.006+02:00

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.

Gestión de plagas en Andalucia: ¿Es efectiva la fumigación de rodales para combatir la procesionaria?

2009-09-28T15:50:00.009+02:00

Anualmente la Consejería de Medio Ambiente de la Junta de Andalucía invierte cerca de un millón o millón y medio de euros en la fumigación de rodales atacados por procesionaria. Sin embargo, cabe preguntarse si dicha medida es efectiva y si la disminución de los brotes de procesionaria responde a una respuesta a la fumigación o, por el contrario, forma parte de los ciclos naturales de esta plaga. Responder a esta pregunta es, por tanto, de gran importancia para la gestión si se quieren minimizar los costes de tratamientos que son, a veces, tan caros como ineficaces.

Para ello, nuestro grupo de investigación ha trabajado con una base de datos recopilada por la Junta de Andalucía (una de las aplicaciones desarrolladas para la Red de Daños y Equilibrios) que contiene información sobre el grado de afectación de los rodales forestales y de los tratamientos aplicados sobre dichos rodales desde 1992 hasta la fecha. Anteriormente se analizó el potencial de una versión preliminar de esta base de datos para el análisis de patrones espacio-temporales de la procesionaria en Andalucía. En este estudio, nos centraremos en la respuesta de la procesionaria en rodales sometidos a fumigación y no tratados. La metodología utilizada es sencilla y consta de los siguientes pasos:

Selección de rodales con un grado de afectación de 3 o más (la escala ordinal utilizada va de 0 a 5, en dónde 3 = defoliación fuerte y bolsones numerosos en los bordes del rodal y algo de defoliación en el centro del rodal) que hayan sido sometidos en el otoño de ese mismo año (septiembre-octubre) a fumigación aérea con objeto de controlar la plaga (los datos más completos se tienen sólo para el período 2002-2005 y es con estos datos con los que se han realizado los sucesivos análisis). Hay que puntualizar que, en principio, sólo se fumigan rodales con un grado de afección de 3 o más, y por eso tomamos este criterio a la hora de seleccionar los rodales para el análisis.
Cálculo de un "índice de recuperación", que se calcularía como la diferencia entre el grado de afectación de un año y el siguiente. Si un año el rodal ha sido asignado a un nivel de daño 3 y en el siguiente año, el grado de afección disminuye, el índice tendría un valor negativo e indicaría una buena recuperación del rodal.
Selección de rodales con un grado de afectación de 3 o más que no hayan sido fumigados. Las muestras se aparean de tal forma que, en cada año, cada rodal con un grado de afectación de 3 o más fumigado se "empareja" con el rodal más próximo geográficamente que también haya tenido un grado de afectación de 3 o más pero que no haya sido fumigado.
Se utiliza un test de la t pareado para comparar los índices de recuperación de rodales tratados y no tratados en cada una de las cuatro especies principales de pino (Pinus halepensis, P. nigra, P. pinaster y P. pinea). Para P. sylvestris no hubo una muestra suficientemente representativa cómo para realizar un test estadístico.

Los resultados mostraron que no hay diferencias significativas en el índice de recuperación de rodales fumigados y no fumigados para las cuatro especies: (Figura 1) Pinus halepensis (p-value = 0.7738), P. nigra (p-value = 0.6987), P. pinaster (p-value = 0.2939) y P. pinea (p-value = 0.4793). Esto indica claramente que, cuando menos, no hay evidencias de que los rodales sometidos a fumigación se recuperen antes que los rodales no fumigados y por tanto cuestiona la validez de esta práctica tan habitualmente usada en la gestión forestal.

Figura 1. Distribuciones del índice de recuperación de rodales forestales fumigados (línea contínua) y no fumigados (línea discontínua). Las barras verticales indican la media de las distribuciones de este índice para cada grupo en cada una de las cuatro especies analizadas: P. halepensis (n = 19 pares de rodales), P. nigra (n = 26 pares de rodales), P. pinaster (n = 14 pares de rodales), P. pinea (n = 71 pares de rodales).

Por otro lado, la Comisión Europea quiere prohibir con carácter general, o como mínimo restringir al máximo, el uso de la fumigación aérea en la UE por los daños que este método puede provocar sobre el medioambiente y la salud de las personas.

Queda claro el mensaje, si no es efectivo y encima causa riesgos innecesarios ¿para qué se siguen fumigando miles de hectáreas cada año?

Acortar las URL largas en documentos

2009-09-16T17:13:00.003+02:00

Un pequeño inconveniente que surge cuando se prepara material didáctico (o informes) es que las URL que se citan (a veces para referirse a datos subidos en una dirección web, o documentos que se pueden descargar on-line) son muy largas. Esto supone un inconveniente si los estudiantes tienen que reescribir toda la dirección, pudiendo dar lugar a errores y perder un tiempo innecesario. Otro problema que he encontrado más recientemente es que en Lyx (ver entrada del 7 de septiembre de 2009), las URL largas quedan fuera de los márgenes del documento, ya que el programa interpreta la URL (o gran parte de ella) como una única palabra. Esto hace que no se lea la dirección completa y que estéticamente quede mal.

Posiblemente esto sea debido a un bug del programa, aunque he tenido exactamente el mismo problema al escribir entradas en el blog. En cualquier caso, aunque estoy seguro de que debe de haber una solución al problema dentro del propio programa, he encontrado una solución alternativa que implica el uso de TinyURL. TinyURL crea una URL pequeña a partir de la URL introducida. Esta URL es mucho más cómoda de manejar y además es permanente. Espero que esta herramienta sea de utilidad a los usuarios. ¡A mi me ha simplificado mucho la vida!

Lyx y Sweave

2009-09-07T17:16:00.005+02:00

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

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

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

Principales diferencias de Lyx con otros editores de texto

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

Cómo instalar Lyx

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

Lyx y Sweave

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

Instalar Sweave

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

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

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

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

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

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

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

2009-07-31T10:25:00.009+02:00

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.

Nuevo artículo publicado en Ecology and Society

2009-07-22T20:26:00.006+02:00

Para todos aquellos que estén interesados en cuestiones relacionadas con el manejo forestal sustentable y los efectos del impacto humano sobre la diversidad forestal, aquí va un artículo muy interesante que acabamos de publicar en Ecology and Society. Este artículo resume los resultados de las investigaciones realizadas por varios grupos europeos y latinoamericanos en el marco de dos proyectos europeos consecutivos del programa INCO (SUCRE y BIOCORES), con especial énfasis en bosques tropicales montanos de Centroamérica y bosques templados de Sudamérica.

El artículo recalca la necesidad de crear herramientas que permitan evaluar de manera integrada el impacto humano sobre la diversidad forestal y que apoyen la toma de decisiones en lo referente al manejo forestal.

Puede descargarse en formato html o pdf.

Diversidad florística y cambios de uso del suelo en Chiapas

2009-07-08T18:59:00.002+02:00

Recientemente ha salido publicado en la revista Investigación Ambiental, Ciencia y Política Pública, un artículo sobre las tendencias y proyecciones del uso del suelo y la diversidad florística en los Altos de Chiapas, México. Este artículo revisa evidencias recientes de la deforestación y de impactos del cambio de uso del suelo en la poblada región conocida como Los Altos o la Meseta Central de Chiapas. Las preguntas guía en torno a las cuales se estructura este trabajo son: (1) ¿qué referencia de riqueza florística es posible proponer para valorar el grado de afectación de los paisajes agropecuarios, que a la vez definiría el potencial y tamaño del reto que representa la restauración ecológica en la región?, (2) ¿cuáles son las principales tendencias de uso del suelo?, (3) ¿cuáles son los principales procesos y consecuencias físicas y biológicas de estos cambios? y (4), ¿cuáles son las tendencias y mayores amenazas para la conservación y aprovechamiento sustentable de los bosques en la región?

Los resultados de este trabajo resaltan la necesidad de definir los elementos de una estrategia para la restauración de áreas deforestadas y el enriquecimiento florístico de los bosques remanentes.

Referencia

González-Espinosa, M., Ramírez-Marcial, N., Galindo-Jaimes, L., Camacho-Cruz, A., Golicher, D., Cayuela, L. y Rey-Benayas, J.M. 2009. Tendencias y proyecciones del uso del suelo y la diversidad florística en los Altos de Chiapas, México. Investigación Ambiental, Ciencia y Política Pública 1(1): 40-53.

Curso de análisis de datos en R - Sesión 5

2009-05-13T13:09:00.003+02:00

Técnicas de análisis de datos multivariantes

En esta clase aprenderemos a manejar distintas técnicas multivariantes en R, como análisis de componentes principales (PCA), análisis de ordenación o análisis de la varianza multivariante (solamente los árboles de regresión y clasificación no serían considerados como una técnica multivariante). El listado completo de técnicas está enumerado a continuación:

Análisis de componentes principales
Análisis de la varianza multivariado (MANOVA)
Escalamiento multidimensional no métrico (NMDS)
Análisis de correspondencias canónico (CCA)
Árboles de regresión y clasificación (CART)

La clase se dividirá en grupos. Cada grupo deberá elegir una técnica determinada, leer la documentación sobre lo que esa técnica hace e imaginar situaciones en sus respectivos campos de investigación en las que el uso de esta técnica podría ser de utilidad. Toda esta información está disponible en una wiki. Una wiki es una herramienta que nos permite trabajar de forma colaborativa on-line.

Una vez leída la documentación, cada grupo deberá aplicar dicha técnica a la resolución de un caso de estudio. Para ello, tendrá que escribir el código necesario que permitirá implementar esa función en R con los datos provistos en el ejemplo. Esto llevará aproximadamente la primera mitad de la clase. En la segunda mitad de la clase, cada grupo explicará a sus compañeros los fundamentos básicos de esa técnica y mostrará su implementación en R.

Curso de análisis de datos en R - Sesión 4

2009-05-11T15:15:00.007+02:00

Modelos Lineales Generalizados (GLM)

Los modelos lineales (regresión, ANOVA, ANCOVA) se basan en los siguientes supuestos:

los errores se distribuyen normalmente;
la varianza es constante; y
la variable respuesta se relaciona linealmente con la(s) variable(s) independiente(s).

En muchas ocasiones, sin embargo, nos encontramos con que uno o varios de estos supuestos no se cumplen. Por ejemplo, es muy común en ecología que a medida que aumenta la media de la muestra, aumente también su varianza. Estos problemas se pueden llegar a solucionar mediante la transformación de la variable respuesta (por ejemplo tomando logaritmos). Sin embargo estas transformaciones no siempre consiguen corregir la falta de normalidad, la heterocedasticidad (varianza no constante) o la no linealidad de nuestros datos. Además resulta muchas veces difícil interpretar los resultados obtenidos. Si decimos que la abundancia de pino silvestre es función de la elevación tenemos una idea más o menos clara de lo que esto puede significar. Si la relación es positiva, un aumento de la elevación aumentaría la abundancia de esta especie. Pero ¿qué quiere decir que el logaritmo de la abundancia de pino silvestre es función de la elevación? Esto ya no es tan intuitivo. La cosa se complica aún más cuando utilizamos otro tipo de transformaciones, como las exponenciales, las potencias, etc. Una alternativa a la transformación de la variable respuesta y a la falta de normalidad es el uso de los modelos lineales generalizados.

Los modelos lineales generalizados (GLM de las siglas en inglés de Generalized Linear Models) son una extensión de los modelos lineales que permiten utilizar distribuciones no normales de los errores (binomiales, Poisson, gamma, etc.) y varianzas no constantes. Ciertos tipos de variables respuesta sufren invariablemente la violación de estos dos supuestos de los modelos normales y los GLM ofrecen una buena alternativa para tratarlos. Específicamente, podemos considerar utilizar GLM cuando la variable respuesta es:

un conteo de casos (p.e. abundancia de una planta);
un conteo de casos expresados como proporciones (p.e. porcentaje de plántulas muertas en un experimento de vivero);
una respuesta binaria (p.e. vivo o muerto, hombre o mujer).

Curso de análisis de datos en R - Sesión 3

2009-05-04T11:40:00.004+02:00

Modelos lineales en R: Regresión, ANOVA y ANCOVA

¿Qué es una regresión? ¿Y un ANOVA? ¿Cuál es la principal diferencia entre ambos? ¿Qué supuestos estadísticos debemos asumir cuando llevemos a cabo este tipo de análisis? Estas y otras preguntas son críticas en la aplicación de modelos lineales a la resolución de problemas estadísticos.

En esta sesión se analizan distintos casos de estudio mediante el uso de modelos lineales y se explica cómo evaluar los supuestos de dichos modelos, cómo solucionar problemas de colinealidad y cómo estandarizar las variables para poder comparar los coeficientes del modelo resultante.

Los pasos a seguir para ajustar un modelo lineal (y prácticamente casi cualquier otro modelo estadístico paramétrico) se resumen en la siguiente figura.

En esta sesión se verán los siguientes contenidos:

Conceptos estadísticos básicos: ANOVA y regresión
Cosas importantes antes de empezar
Cómo ajustar un modelo lineal en R

Un ejemplo de regresión
Un ejemplo de ANOVA
Un ejemplo de ANCOVA
Interacción entre factores o factores y co-variables

Evaluación de los supuestos del modelo: Exploración de los residuos
Problemas de colinealidad: Reducción de variables
Estandarización de coeficientes

Curso de análisis de datos en R - Sesión 2

2009-05-01T11:25:00.003+02:00

Gráficos en R

R ofrece una increíble variedad de gráficos. Para tener una idea, escribe el comando demo(graphics). Cada función gráfica en R tiene un enorme número de opciones permitiendo una gran flexibilidad en la producción de gráficos y el uso de cualquier otro paquete gráfico palidece en comparación. Al contrario que con funciones estadísticas, el resultado de una función gráfica no puede ser asignado a un objeto sino que es enviado a un dispositivo gráfico. Un dispositivo gráfico es una ventana gráfica o un archivo.

Existen dos tipos fundamentales de funciones gráficas: las funciones gráficas de alto nivel que crean una nueva gráfica y las funciones gráficas de bajo nivel que agregan elementos a una gráfica ya existente. Las gráficas se producen con respecto a parámetros gráficos que están definidos por defecto y pueden ser modificados con la función par().

A lo largo de esta sesión el alumno aprenderá a manejar gráficos y dispositivos gráficos en R. También se verá en detalle algunas funciones gráficas y sus parámetros, ejemplos prácticos del uso de estas funciones en la producción de gráficos y la descripción de los distintos paquetes que las contienen. Los contenidos concretos se detallan a continuación:

Conceptos básicos
La organización de R

Tipos de sistemas gráficos
Formato de gráficos

Gráficos tradicionales

El paquete graphics
Funciones de alto nivel: representación de una o dos variables
Funciones de alto nivel: representación de múltiples variables
Funciones de bajo nivel
Cómo personalizar un gráfico

Los paquetes grid y lattice

Curso de análisis de datos en R - Sesión 1

2009-04-26T10:21:00.006+02:00

Una introducción a R

R es un lenguaje orientado a objetos. Aunque existen algunas interfaces gráficas para R como Rcommander, es muy recomendable aprender R como un lenguaje en vez de tratarlo como un programa estadístico convencional. Como ambiente de trabajo, R ofrece una serie de ventajas:

Sus posibilidades gráficas son excelentes
Es muy flexible. Los procedimientos estadísticos estándar se pueden aplicar con sólo utilizar el comando apropiado. Además, existen multitud de librerías (a los que llamaremos paquetes de ahora en adelante) programadas por los usuarios de todo el mundo para llevar a cabo procedimientos específicos.
Es libre. Libre no quiere decir gratuito (aunque R también lo es). Libre significa que podemos acceder al código escrito por otros usuarios y modificarlo libremente. A pesar de que R viene sin garantía alguna (al iniciar la sesión de R saldrá la siguiente advertencia “R is free software and comes with ABSOLUTELY NO WARRANTY”), la mayor parte del código de R, o por lo menos, el código más comúnmente utilizado por los usuarios, ha sido meticulosamente supervisado por estadísticos y académicos de mucho prestigio de todo el mundo (el llamado “R Core team”).
Podemos además programar nuestras propios procedimientos y aplicaciones
En la misma página desde la que se puede bajar el programa, existe abundante documentación sobre cómo utilizarlo
Es gratuito.

Esta documentación forma parte del Curso de Doctorado de Análisis de Datos Ecológicos en R, que se impartirá entre el 5 y el 13 de mayo de 2009 en el Departamento de Ecología de la Universidad de Alcalá. La primera sesión constituye una introducción somera a R, en dónde el alumno podrá explorar los conceptos más básicos del manejo de esta herramienta, desde la instalación del software a la manipulación de datos y la programación de funciones sencillas. La documentación para esta primera sesión se puede descargar aquí. El índice de contenidos se detalla a continuación:

¿ Qué es R?
¿Cómo instalar R?
CRAN y paquetes
Tipos de objetos en R y la función str()
Funciones y argumentos
El menú de ayuda: Aprendiendo a ser autosuficientes
Como leer datos en R
Recomendaciones para organizar una sesión de trabajo
Funciones básicas para la manipulación de datos
Como repetir un procedimiento con el comando for()
Citando R en los trabajos científicos
Referencias

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

2009-04-13T14:32:00.008+02:00

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:

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.
Se calculó la integral total anual (NDVI-I) como la suma del área bajo la curva en cada intervalo.
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.

Una nueva revista de conservación tropical: Tropical Conservation Science

2009-04-13T13:38:00.004+02:00

Hace poco descubrí una nueva revista de investigación centrada en la conservación de los ecosistemas tropicales: Tropical Conservation Science. Esta revista, que se viene publicando desde el año 2008 (tan sólo cuenta con cinco volúmenes publicados hasta el momento), contiene artículos muy interesantes sobre distintos aspectos relacionados con la conservación de los ecosistemas tropicales, como la política territorial, la conservación de comunidades y especies, o la conservación genética de poblaciones, entre otros muchos.

Tropical Conservation Science es una revista de libre acceso (open source), lo que favorece la diseminación de los trabajos científicos frente a revistas de suscripción. El Dr. Alejandro Estrada, un eminente primatólogo de la Universidad Nacional Autónoma de México, es junto con Rhett Butler, de mongabay.com, uno de los dos editores ejecutivos de la revista. La revista cuenta además con una larga lista de editores asociados de la talla y calidad del Dr. Rodolfo Dirzo (Stanford University, USA), la Dra. Sarah Boyle (Arizona State University, USA) o el Dr. William F. Laurance (The Smithsonian Tropical Research Institute, Panama), entre otros.

Tropical Conservation Science ofrece, por tanto, una alternativa sólida y con futuro de proyección a revistas científicas ya consolidadas en el ámbito de la ecología tropical, como Biotropica o Journal of Tropical Ecology. Por todo ello, espero que esta breve reseña contribuya a su difusión, para que sea conocida por un público amplio y llegue a ser en el futuro una revista de referencia en el campo de la conservación de los ecosistemas tropicales.

Presentación al 5º Congreso Forestal Español

2009-03-30T11:15:00.014+02:00

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

2009-03-17T14:20:00.016+01:00

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:

Porcentaje de decaimiento observado in situ en campo.
Especie (Pinus nigra y P. sylvestris).
Área basimétrica del rodal (var. estructural).
Densidad de pies del rodal (var. estructural).
Coordenadas geográficas x e y (var. espaciales).
Elevación (var. topográfica).
Pendiente (var. topográfica).
Orientación (var. topográfica).
Insolación en el solsticio de invierno (var. topográfica).
Insolación en el solsticio de verano (var. topográfica).
Í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.

Como utilizar el potencial de GRASS desde R

2009-03-09T17:34:00.011+01:00

R tiene varios paquetes que pueden utilizarse para el análisis y visualización de datos espaciales (sp, maptools, rgdal, etc.). Sin embargo, R no es un SIG y, por tanto, tiene ciertas limitaciones a la hora de trabajar con información espacial. Por ejemplo, para algo tan simple como extraer la pendiente y la orientación a partir de un modelo digital de elevaciones, R necesita del paquete RSAGA, que a su vez necesita del software gratuito SAGA diseñado por Alexander Brenning, que funciona únicamente bajo Windows. Esto supone ciertas limitaciones al uso de R como SIG. Estas limitaciones no existen cuando se utiliza GRASS en combinación con R. GRASS ofrece un complemento perfecto a R para la implementación de operaciones típicamente realizadas en SIG y R permite potenciar las capacidades de GRASS para la implementación de análisis geoespaciales.

Para poder hacer uso de ambos programas en necesario tener instalados GRASS y una versión de R posterior a 2.1.0. Desde el 'shell' se puede abrir una sesión de R. Lo primero será instalar los paquetes necesarios para poder utilizar R en combinación con GRASS (sp, rgdal, maptools, spgrass6, spGDAL, spmaptools).

GRASS 6.3.0 (Filabres):~ > R

install.packages(c("sp", "rgdal", "maptools"), dependencies= TRUE)

rS <- "http://r-spatial.sourceforge.net/R" install.packages(c("spgrass6", "spGDAL", "spmaptools"), repos=rS, dependencies=TRUE)

Sin salir de la sesión de R, debemos de cargar el paquete spgrass6 y la 'location' de GRASS donde se encuentran todas las capas de información raster o vectorial.

library(spgrass6)
G <- gmeta6()
str(G)

Si ahora queremos acceder a alguna de las capas raster del 'location' en el que estamos trabajando dentro de GRASS tendremos que utilizar la función de R readCELL6sp() ò readRAST6(). Esto genera un objeto del tipo SpatialGridDataFrame.

mde <- readCELL6sp("mde")

Sin embargo, es importante tener en cuenta que si cargamos muchas capas simultáneamente, especialmente si son capas con un gran número de celdas, se puede agotar la memoria virtual. En este caso R dará el mensaje 'vector memory exhausted'. Cuando esto me ha ocurrido a mi, R se ha quedado colgado y he tenido que reiniciar la sesión. Así que lo que aconsejo es ir eliminando los objetos que ya no necesitemos una vez que extraigamos de ellos la información necesaria. Para conocer el uso que se va haciendo de la memoria se puede utilizar la función gc().

gc()

used (Mb) gc trigger (Mb) max used (Mb)
Ncells 279590 7.5 531268 14.2 361111 9.7
Vcells 3160949 24.2 15608917 119.1 30379099 231.8

Es decir que se están utilizando alrededor de 31 MB de memoria virtual (hay que mirar la primera columna). Ahora eliminamos la capa mde.

rm(mde)
gc()

used (Mb) gc trigger (Mb) max used (Mb)
Ncells 279513 7.5 531268 14.2 361111 9.7
Vcells 137524 1.1 12487133 95.3 30379099 231.8

Y comprobamos que disminuye considerablemente el uso de la memoria virtual.

Antes de acceder a las capas raster, vamos a leer una capa de puntos en formato *.txt (parcelas muestreadas en campo) para la cuál extraeremos la información de elevación, pendiente, altitud e insolación de las capas raster.

points <- read.table(url("http://archivos-para-subir.googlegroups.com/web/Analisis+correlacionales.txt?gda=0ywgME4AAACkOfwqHAVd4YqgfIB09GDRN3No94DJejqv6LWrffOcUtQTec6xGbbz9B_IWZmp37itXp5Ud9d8afGj09bAbTSQ47Cl1bPl-23V2XOW7kn5sQ"), header = T, sep = "\t", dec = ",")
str(points)

'data.frame': 76 obs. of 7 variables:
$ x : num 543726 543486 542903 535211 539086 ...
$ y : num 4124164 4124224 4122417 4118071 4121908 ...
$ Id : Factor w/ 76 levels "1_12","1_13",..: 1 2 4 6 7 8 9 10 11 12 ...
$ Especie : Factor w/ 2 levels "Pinus nigra",..: 1 1 1 1 1 1 1 2 1 1 ...
$ Defoliacion: num 55 55 30 55 30 55 55 60 45 55 ...
$ Area_basim : num 25.8 28.1 24.2 22.0 18.0 ...
$ Densidad : num 286 286 1592 923 1974 ...

El siguiente paso es convertir esta capa de puntos en un objeto espacial del tipo SpatialPointsDataFrame y asignarle una proyección (ha de ser la misma que la de las capas raster y sino habría que reproyectar esta capa).

coordinates(points) <- data.frame(points$x, points$y)
proj4string(points) <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")

Y ahora sí, importamos cada una de las capas raster de GRASS a R, extraemos los valores de las capas raster para los puntos incluidos en la capa vectorial points y, por último, borramos las capas raster para no agotar la memoria virtual.

mde <- readCELL6sp("mde")
extract.mde <- overlay(mde, points)
rm(mde)

slope <- readCELL6sp("slope")
extract.slope <- overlay(slope, points)
rm(slope)

aspect <- readCELL6sp("aspect")
extract.aspect <- overlay(aspect, points)
rm(aspect)

rad.summer <- readCELL6sp("beam_radiation_172")
extract.rad.sum <- overlay(rad.summer, points)
rm(rad.summer)

rad.winter <- readCELL6sp("beam_radiation_354")

extract.rad.win <- overlay(rad.winter, points)
rm(rad.winter)

wetness <- readCELL6sp("wetness_index")
extract.wetness <- overlay(wetness, points)
rm(wetness)

Y ahora juntamos toda la información extraida en un único data.frame.

points@data <- cbind(points@data, mde = extract.mde@data, slope = extract.slope@data, aspect = extract.aspect@data, rad.summer = extract.rad.sum@data, rad.winter = extract.rad.win@data, wetness = extract.wetness)

El resultado final, es por tanto, un data.frame con información sobre valores de decaimiento (incluidos en el archivo points) y toda una serie de variables físicas obtenidas a partir de capas raster (pendiente, altitud, etc) que podemos analizar estadísticamente en R, cosa que en principio no podríamos hacer en GRASS. En la siguiente entrada se muestra un análisis preliminar de estos datos en R.

Cálculo de superficies de insolación

2009-03-06T17:10:00.011+01:00

En GRASS es relativamente sencillo implementar un modelo para estimar la insolación a partir de los mapas de altitud, pendiente y orientación por medio del comando r.sun.

Se puede acceder a este comando a través del 'shell' o del gestor GIS. Si optamos por la segunda opción, hay que ir a Raster - Solar radiance and shadows - Solar irradiance and irradiation. Se abrirá una ventana con tres pestañas: (1) Options; (2) Input_options; (3) Output_options. Algunos de los campos que hay que rellenar son obligatorios pero muchos otros son optativos y, o bien no son completamente necesarios, o están definidos por defecto.

Antes de utilizar este comando es conveniente leerse la ayuda de la función. r.sun puede utilizarse para calcular tanto la radiación directa, difusa y reflejada [W.m-2] en un momento puntual del día (modo 1), como la radiación directa, difusa y reflejada acumulada [Wh.m-2.day-1] a lo largo de todo un día (modo 2). El segundo modo suele ser más informativo, sobretodo en estudios relacionados con la producción de cultivos o la presencia de formaciones vegetales. Esta radiación acumulada será la que calculemos en este ejemplo.

En la primera pestaña (Options) se pide información sobre:

el día del año en el que se quiere calcular la radiación (día juliano).
el intervalo del tiempo en el que se calcula la radiación cuando se estima la radiación acumulada durante todo el día. Esto es sólo aplicable al modo 2. 0.5 (valor por defecto) significaría que se estima la radiación cada media hora a lo largo del día y luego se calcula la integral para todas las estimas puntuales.
el valor de declinación solar (radianes) se estima por defecto de manera automática utilizando el día del año. No obstante, si se tiene una medida más precisa de este valor, puede introducirse en la casilla correspondiente.
la hora solar para la que se quiere calcular la radiación, sólo aplicable en el caso del modo 1.

Por tanto, en el caso de estar interesados en el cálculo de la radiación acumulada (modo 2), sólo se tiene que especificar el día del año y el intervalo de tiempo para el cuál se hacen las estimas parciales.

En la siguiente pestaña (Input_options) hay que introducir el nombre de los archivos que contienen las capas raster (en formato Grass) de elevación, (elevation) pendiente (slope) y orientación (aspect).

Por defecto, hay establecido un coeficiente de turbidez atmosférica (Linke atmospheric turbidity coefficient) de 3.0. Este coeficiente refleja la presencia de aerosoles y nubosidad que puede influir en la radiación difusa y reflejada. En la página de ayuda muestran una tabla en donde se sugieren coeficientes de Linke orientativos según el tipo de área de estudio y mes del año. Si se dispusiera de una capa raster con el valor de este coeficiente para toda el área de estudio, se podría introducir en la casilla que dice 'Name of the Linke atmospheric turbidity coefficient input raster map'. Como generalmente no se dispone de una información tan detallada, se deja este espacio en blanco.

En la tercera pestaña (Output_options) es necesario especificar el nombre de los archivos de salida para la radiación directa (beam_rad), difusa (diff_rad) y reflejada (refl_rad). Además, en el caso del modo 1 habrá que especificar el nombre de un archivo de salida con el ángulo de incidencia del sol (incidout). En el caso del modo 2, habrá que especificar el nombre de un archivo de salida con el tiempo de insolación (insol_time), es decir, el número de horas a lo largo del día que cada pixel ha estado expuesto a la radiación directa.

Una vez definidos todos estos parámetros, le damos a RUN. El proceso, especialmente en el caso del modo 2, puede ser muy lento ya que exije una gran cantidad de cálculos computacionales. Para el caso de la Sierra de los Filabres (Almería), el cálculo de la insolación diaria total (radiación integrada) para un área de 3839 x 7086 píxeles de 10 x 10 m duró aproximadamente 3-4 horas.

Por otro lado hay que tener en cuenta que para la definición de las características de insolación de un lugar deben estimarse los modelos de insolación representativos de un conjunto suficiente de periodos anuales. En efecto, al menos en latitudes medias o altas los contrastes estacionales son fuertes y contienen una información ambiental que no puede despreciarse. Por este motivo no es suficiente caracterizar la insolación para un valor medio anual, al menos en zonas no ecuatoriales. Por esta razón, he calculado la insolación en el solsticio de verano (día juliano 172) y en el solsticio de invierno (día juliano 354). Luego estas dos variables pueden considerarse individual o conjuntamente (mediante la media) en análisis posteriores.

Por último, es posible obtener una representación del mapa de insolación utilizando el siguiente código en la consola de comandos de GRASS.

g.region rast=mde
d.mon x1
r.colors beam_radiation_172 col=grey
d.his i_map=beam_radiation_172 h_map=mde

Análisis de modelos digitales de elevación en GRASS

2009-03-06T12:43:00.007+01:00

En GRASS se pueden calcular los mapas de pendiente y orientación a partir de un Modelo Digital de Elevaciones (Digital Elevation Model (DEM) en inglés) desde el administrador GIS pulsando:

Raster - Terrain analysis - Slope and aspect

O bien desde la consola de comandos (shell) con el comando r.slope.aspect.

r.slope.aspect el=dem as=aspect sl=slope

Los dos mapas son procesados a la vez. Los ángulos del mapa de orientación son medidos en el sentido opuesto a las manecillas del reloj a partir del Este. Las pendientes son calculadas por defecto en grados.

Si se quieren visualizar estos dos archivos basta con usar los comandos d.mon y d.rast para abrir un monitor y visualizar (display) en él la capa raster respectivamente. Alternativamente el comando d.rast.leg agrega una leyenda sencilla al monitor.

d.mon x1
d.rast aspect
d.mon x2
d.rast.leg slope

Existen módulos adicionales para trabajar con los DEMs: las áreas con depresiones pueden ser llenadas con r.fill.dir, las líneas de flujo pueden ser calculadas con r.flow. El análisis de cuencas puede ser realizado con r.watershed y con r.terraflow en archivos muy grandes. Todos estos comandos pueden ser utilizados también desde el administrador GIS (en la pestaña raster).

¡Bienvenido a GRASS!

2009-03-06T09:53:00.024+01:00

GRASS (por su siglas en inglés, the Geographic Resources Analysis and Support System, http://grass.itc.it) es un Sistema de Información Geografica (SIG) gratuito y de código abierto, integrado con un sistema de procesamiento de imágenes y subsistemas de visualización. GRASS permite trabajar con imágenes raster y vectorial. Tiene interfaces paraPostgreSQL, MySQL, DBF y bases de datos conectadas con ODBC. Además puede ser conectado a UMN/MapServer, R, gstat, Matlab, Octave, Povray y otros software.

Hasta hace poco siempre había trabajado con otros software comerciales como ArcGIS, Idrisi, PCI o ERDAS. GRASS siempre me había parecido muy complicado y poco amigable. Sin embargo, sus características como software libre y gratuito lo convierten en una herramienta muy versátil y potente, por lo que merece la pena hacer el esfuerzo de aprendizaje. En esta entrada explico como empezar a trabajar en GRASS y los conceptos más básicos para aquellos que, como yo hasta ayer, quieran dar el salto pero no acaban de animarse.

Instalación y manuales

GRASS está disponible desde el ITC-irst en Italia (http://grass.itc.it). GRASS está diseñado para correr en varios ambientes UNIX, tales como GNU/Linux, SUN-Solaris, Irix o MacOS X. También puede correr bajo MSWindows/NT/2000/XP con Cygwin (un simulador de Linux), aunque posiblemente las próximas versiones de GRASS se puedan correr bajo Windows directamente. Libros, tutoriales, manuales, cursos en línea y mas documentación se encuentra listada en el 'Proyecto para la Documentación de GRASS' (http://grass.itc.it/gdp). En español, hay al menos un par de tutoriales de introducción a GRASS (http://mpa.itc.it/markus/osg05/). Ambos son traducciones de un tutorial de Markus Neteler para un taller de trabajo y no difieren apenas en cuanto a contenidos entre sí. En mi opininión, estos tutoriales no son particularmente buenos, pero por lo menos sirven para entender algunos conceptos básicos del funcionamiento de GRASS, por lo que merece la pena echarles un vistazo (muchos de los contenidos de esta entrada han sido extraidos de estos tutoriales).

Estructura de los proyectos en GRASS

Cuando alguien está familiarizado con el uso de SIG y sensores remotos pero nunca antes había trabajado con GRASS, lo más importante es saber cómo se estructuran los proyectos y cómo abrir una sesión en GRASS. Esto es lo que normalmente más cuesta entender. La información de GRASS se encuentra almacenada en un directorio conocido como 'database' (también llamado 'GISDBASE') . Este directorio se debe de crear antes de empezar a trabajar con GRASS (yo lo he llamado 'home/luis/grassdata'). Dentro de este 'database', los proyectos están organizados por áreas de proyectos en subdirectorios llamados 'locations'. Un 'location' está definido por un sistema de coordenadas, una proyección cartográfica y unas fronteras geográficas. Los subdirectorios y archivos que definen un 'location' son creados automáticamente cuando GRASS es iniciado la primera vez con un nuevo 'location' (esto es lo que vamos a hacer un poco más abajo). Cada 'location' puede tener múltiples 'mapsets'. Un motivo para mantener múltiples 'mapsets' es almacenar mapas relacionados a algún aspecto específico del proyecto o subregiones. Otro motivo es permitir el trabajo simultáneo de usuarios a capas almacenadas dentro del el mismo 'location', por ejemplo, equipos trabajando en elmismo proyecto. Para trabajos en equipo, una base de datos de GRASS centralizada puede ser creada en un sistema de archivos de red (por ejemplo, NFS). Además de acceder el propio 'mapset', cada usuario puede leer capas de otros 'mapsets', aunque cada usuario solo puede modificar o borrar las capas en su propio 'mapset'.

Cuando se crea un nuevo 'location', GRASS automáticamente crea un 'mapset' especial llamado PERMANENT donde la información principal del proyecto puede ser almacenada. Datos en el 'mapset' PERMANENT pueden ser solamente agregados, modificados o borrados por el propietario del 'mapset' PERMANENT. Sin embargo, los datos pueden ser consultados, analizados y copiados por otros usuarios a sus respectivos 'mapsets'. El 'mapset' PERMANENT es útil para entregar información general (ejemplo, un modelo de elevación digital), accesible pero protegida contra escritura de todos los usuarios que se encuentran trabajando en el mismo 'location' del propietario de la base de datos. En principio si no se tiene pensado trabajar en equipos con varios usuarios y con una gran cantidad de información, se puede trabajar directamente en el 'mapset' PERMANENT.

Primeros pasos en GRASS: Cómo crear un nuevo 'location'

Una gran diferencia entre GRASS y otros SIG es que GRASS requiere los parámetros de la proyección antes de que el usuario pueda trabajar en un 'location'. La ventaja de esto es que las cosas quedan bien definidas y se evita un desorden con la combinación de proyecciones, cosa que en mi experiencia es bastante común cuando se trabaja con otros software que no siguen esta estructura de proyectos. Para crear un nuevo 'location' primero tenemos que abrir GRASS. En Linux (yo trabajo actualmente con Ubuntu Intrepid Ibex 8.10) hay que ir a la terminal (konsole) y escribir:

grass63

o simplemente:

grass

Esto abrira la ventana de bienvenida de GRASS.

Lo primero que hay que hacer es seleccionar el 'database' donde se va a guardar el nuevo 'location'. Este directorio se ha tenido que crear antes de empezar a trabajar en GRASS. Ahora selecciona el botón 'Projection values' para crear un nuevo 'location' (Define new location with...), el cual te llevará a una nueva ventana de texto.

En la pantalla ingresa un nuevo nombre para el 'location' (no puede contener espacios en blanco) y luego continua presionando la tecla ESC y la tecla RETURN de tu teclado. Te aparecerá otra ventana de texto diciendo que esa 'location' no existe, cuáles son las 'locations' disponibles y preguntándote si quieres crear una nueva 'location' con el nombre asignado. A esto último se dice que sí y se pulsa RETURN.

Esto te lleva a otra ventana de texto en dónde te explican la información que vas a necesitar para crear un nuevo 'location' (sistema de coordenadas, zona UTM y parámetros de la proyección, coordenadas del área de trabajo, etc). Te preguntará si tienes toda esa información. Se dice que sí y se pulsa RETURN.

Una nueva ventana te pedirá que definas el sistema de coordenadas para tu 'location'. En mi caso defino C (UTM) y pulso RETURN.

Lo siguiente es ingresar una pequeña descripción o título para el 'location'. En esta misma ventana te pregunta si quieres especificar un datum. Por defecto se responde que sí. Se pulsa RETURN para continuar.

Para especificar el datum se puede pulsar 'list' y te viene una lista de todos los datum disponibles. Para salir de este listado se pulsa 'q'. Se escribe el nombre del datum y se pulsa RETURN.

El último paso es definir las coordenadas del área de estudio. En mi caso por ejemplo, las hice coincidir con las del modelo digital de elevaciones con el que voy a trabajar. En otros casos pueden ser más o menos amplias y abarcar toda una región geográfica como Sierra Nevada, Andalucía o la Península Ibérica. Hay que tener en cuenta que esta es la región que define el 'location', pero que luego es posible redefinir el área de estudio en distintos 'mapset' dentro de un 'location'. Igualmente se pueden redefinir los parámetros de un 'location' una vez que ha sido creado.

Finalmente se pide que se defina la resolución de los píxeles. Una cosa buena de GRASS es que se puede cambiar el tamaño de los píxeles con relativa facilidad. Al contrario que con otros software, no se remuetrea una capa raster, sino que se reproyecta todo un proyecto ('location') a la resolución deseada. Esto agiliza mucho el procesamiento de la información sin necesidad de generar infinidad de nuevos archivos con información parcialmente redundante. Se pulsa RETURN para acabar.

Ya tenemos el proyecto creado. Ahora podemos comenzar GRASS desde la página de bienvenida (paso 1) seleccionando el 'database' donde has creado el nuevo 'location'. (/home/luis/grassdata), el 'location' que acabas de crear (Filabres) y el 'mapset' (PERMANENT). Pulsa 'Enter GRASS'.

¡Y ya podemos empezar a trabajar en GRASS!

Interfaz gráfica y consola de comandos ('shell')

En GRASS se puede trabajar de dos maneras (alternativa o simultáneamente): a través de la interfaz gráfica o a través de la consola de comandos ('shell'). Mientras que la primera es más intuitiva, la segunda da mayor versatilidad y potencia a la hora de programar funciones y procesos. Ambas opciones están disponibles todo el tiempo.

La interfaz gráfica consta de tres ventanas.

La ventana principal es el administrador GIS (GIS Manager, izquierda). Desde esta ventana se pueden realizar la mayoría de las funciones relacionadas con el procesamiento de capas raster y vectorial, incluyendo la importación y exportación de capas en otros formatos. También se pueden modificar las características de configuración del 'location'.

Una segunda ventana (Output - GIS.m, arriba derecha) nos muestra los resultados de las operaciones realizadas a través del administrador GIS así como el código implementado (este código se puede usar directamente en la consola de comandos).

Finalmente hay una tercera ventana (Map Display X, abajo derecha) que muestra los mapas seleccionados en el administrador GIS. Otra de las ventajas de GRASS es que en vez de desplegar los mapas en la pantalla, envía el contenido del monitor a un archivo *.PNG por lo que ocupa menos memoria y es mucho más rápido a la hora de visualizar datos, especialmente cuando se trata de capas raster con una gran cantidad de información.

La consola de comandos ('shell') funciona de manera similar a como funcionan la mayoría de las aplicaciones en sistemas Unix. A través del 'shell' podemos obtener toda la funcionalidad que se obtiene a través de la interfaz gráfica de GRASS, pero además podemos programar scripts y conectar con otros software como R.

Importar datos a GRASS

Ya hemos creado nuestra 'location'. Ahora vamos a importar una capa raster a nuestro 'location' para trabajar con ella desde GRASS. GRASS es compatible con la mayoría de los formatos típicos utilizados en SIG y teledetección. En este caso importaremos una capa raster ascii de ArcGIS (GRID). Para ello vamos a:

File - Import raster map . ESRI grid

Alternativamente desde la consola de comandos podríamos haber escrito el comando:

r.in.arc

Cualquiera de las dos formas nos llevarían a la siguiente ventana, en dónde debemos de definir el archivo de entrada y el nombre del archivo de salida que se guardará directamente en nuestro 'location' y 'mapset' definidos. Una vez definidos los argumentos necesarios para ejecutar la función, pulsamos RUN.

Como podemos observar en la parte inferior de la ventana aparece el código que ejecuta todo el proceso incluyendo la definición de los archivos de entrada y salida. Esta línea de texto puede ser ejecutada en la consola de comandos para realizar esta misma operación sin necesidad de tener que abrir un GUI (Graphic's User Interface).

Ahora ya podemos trabajar con esta capa raster en GRASS para realizar distintas operaciones. Para más información sobre cómo cargar capas en GRASS y visualizarlas se puede consultar el manual de Markus Neteler traducido por Guillermo Martínez.

Modelos binomiales en R

2009-01-26T16:28:00.014+01:00

Muchos problemas estadísticos implican el análisis de una respuesta binaria. Por ejemplo, muchas veces se clasifican los individuos muestrales como:

Vivos o muertos,
Presencia o ausencia de especies,
Hombre o mujer,
Sano o enfermo,
Deforestado o no deforestado.

En estos y otros casos, los modelos binarios suponen una herramienta de análisis muy interesante. Para ello es conveniente que, al menos, una de las variables explicativas sea continua. En el caso de la deforestación (o cualquier otro cambio de uso binario), los datos consistirán, por ejemplo, en una lista de variables (ej. distancia a caminos, elevación, pendiente) medidas en píxeles que han sido deforestados en un periodo determinado, y una lista similar para píxeles que no han sido deforestados en ese mismo periodo. La cuestión es tratar de averiguar cuál de estas variables, si es que hay alguna, aumentan la probabilidad de que un píxel individual sea deforestado para ese periodo.

La variable respuesta contiene sólo 0s y 1s. Por ejemplo, 0 representaría píxeles que permanecen como bosque y 1 píxeles que han sido deforestados. Por tanto, habría una única columna para la variable respuesta. Los pasos a seguir para modelar una variable respuesta de tipo binomial son los siguientes:

Crear un vector de 0s y 1s como variable respuesta.
Ajustar un modelo lineal generalizado (funcion glm) con una familia de distribución de errores de tipo binomial (argumento family = binomial).
Reducir el modelo eliminando los términos que no sean significativos y comparando el cambio en devianza con respecto al modelo total.
Interpretación de los resultados.

La función de vínculo en los glm hace referencia a la transformación de la(s) variable(s) antes de ajustar el modelo. En el caso de los modelos binomiales, la función de vínculo que se usa habitualmente es la logística (logit), que transforma con logaritmos naturales la variable respuesta. Por defecto, no hace falta especificar ésto en los modelos.

A la hora de comparar modelos se utiliza la devianza como un indicador de la bonanza del modelo. Para aquellos que no estén muy familiarizados con la jerga estadística, la devianza es una medida del ajuste del modelo en relación al número de parámetros del mismo. Si un modelo explica lo mismo que otro, pero es más sencillo (menor número de párametros) es, por regla general, mejor modelo. Y esto es lo que va a indicar la devianza.

Un ejemplo: El caso de la deforestación en Chile Central

En este ejemplo, la variable respuesta va a ser la ocurrencia de deforestación en Chile Central para el periodo 1975-1985, en dónde 1 = deforestación y 0 = no deforestación. Para obtener estos datos, se llevó a cabo un muestreo sistemático cada 1 km en el área de estudio y se registró la cobertura en 1975 y 1985 para todos aquellos píxeles que eran bosque en 1975. Se midieron también toda una serie de variables explicativas que se detallan a continuación:

Elevación (m),
Pendiente (grados),
Orientación (grados),
Distancia a ríos (m),
Distancia a ciudades grandes (m),
Distancia a pueblos (m),
Distancia a carreteras principales (m),
Distancia a carreteras secundarias (m),
Distancia a suelo agrícola (m),
Distancia al borde del bosque (m).

Los datos originales están en Excel. Para leerlos directamente en R suele ser más rápido pasarlos a *.txt, aunque R puede leer también datos en Excel. Los datos están accesibles en una página web, por lo que se pueden leer desde R utilizando el comando url(). Basta con copiar el código que se ve en azul y pegarlo en R. Funcionará.

fnf75_85 <- read.table(url("http://archivos-para-subir.googlegroups.com/web/F-NF1975-85.txt?gda=f500bkEAAACkOfwqHAVd4YqgfIB09GDRCcOUO-z8uvHRarwVOPM5xOEDwyLfg5d1mmtcx5AUJ3JTCT_pCLcFTwcI3Sro5jAzlXFeCn-cdYleF-vtiGpWAA"), header = T, dec = ".", sep = "\t")

Ya se ha creado un objeto en R llamado fnf75_85, que contiene los datos de las distintas variables que se utilizarán en el análisis. Con el argumento header = T se especifica que la primera línea contiene los nombres de las variables. Con el argumento dec = "." se dice cómo vienen especificados los decimales (alternativamente podría ser ","). Por último, el argumento sep = "\t" especifica que los campos están separados por tabulaciones. Con la función str() se obtiene la estructura del objeto fnf75_85 e información sobre la naturaleza (integer, factor, character, etc) de cada una de las variables.

str(fnf75_85)

'data.frame': 1980 obs. of 15 variables:
$ ID_ALL : int 10 20 77 80 81 118 158 208 255 262 ...
$ ID_FOREST : int 1 2 3 4 5 6 7 8 9 10 ...
$ Class1975 : int 1 1 1 1 1 1 1 1 1 1 ...
$ Class1985 : Factor w/ 8 levels "1","2","3","5",..: 4 3 2 2 2 3 1 1 2 1 ...
$ Change : Factor w/ 3 levels "0","1","NoData": 2 2 2 2 2 2 1 1 2 1 ...
$ Elevation : num 435 136 299 300 300 ...
$ Slope : num 19.6 0.4 0.2 0 0 0 0.2 12.3 0.3 3.8 ...
$ Aspect : num 149 290 357 -1 -1 ...
$ Dist_river : num 4536 1025 751 2005 2672 ...
$ Dist_City.20T: num 27750 52254 28529 25994 25180 ...
$ Dist_village : num 2271 2620 1142 1084 742 ...
$ Dist_road_P : num 4601 2313 4621 1973 1442 ...
$ Dist_road_S : num 752 395 323 552 0 ...
$ Dist_agri : num 94.9 60.0 30.0 30.0 30.0 ...
$ Dist_f_forest: num 84.9 524.8 67.1 161.6 67.1 ...

La variable respuesta es $Change. Se observa que es un factor con 3 niveles: "0", "1" y "NoData". No interesa tener muestras con "NoData", por lo que se deben de eliminar estos datos. Se puede hacer sobreescribiendo el archivo fnf75_85 (! indica exclusión de casos).

fnf75_85 <- fnf75_85[!fnf75_85$Change == "NoData", ]
fnf75_85$Change <- as.integer(fnf75_85$Change)
fnf75_85$Change <- ifelse(fnf75_85$Change == 1, 0, 1)

Y ahora un resumen de la variable respuesta.

table(fnf75_85$Change)

0 1
774 1184

Así que hay 774 píxeles que continúan siendo bosque en 1985 y 1184 píxeles que han sido deforestados entre 1975 y 1985. El diseño es claramente no balanceado (distinto número de casos en la variable respuesta), aunque esto no tiene por qué influir necesariamente en el ajuste del modelo. El ajuste del modelo completo (excluyendo la variable Dist_f_forest, que ha sido mal calculada) se haría de la siguiente forma:

model.full <- glm(Change ~ Elevation + Slope +Aspect + Dist_river + Dist_City.20T + Dist_village + Dist_road_P + Dist_road_S + Dist_agri, data = fnf75_85, family = binomial)

Para explorar los resultados se puede usar la función summary().

summary(model.full)

Call:
glm(formula = Change ~ Elevation + Slope + Aspect + Dist_river +
Dist_City.20T + Dist_village + Dist_road_P + Dist_road_S +
Dist_agri, family = binomial, data = fnf75_85)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.5265 -1.3419 0.9214 1.0039 1.3176

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.177e-01 1.650e-01 3.742 0.000182 ***
Elevation 3.529e-04 2.209e-04 1.597 0.110201
Slope 3.913e-03 4.902e-03 0.798 0.424748
Aspect -8.098e-04 4.641e-04 -1.745 0.081044 .
Dist_river 1.350e-05 6.302e-05 0.214 0.830348
Dist_City.20T -3.136e-06 4.025e-06 -0.779 0.435960
Dist_village 2.353e-05 3.188e-05 0.738 0.460615
Dist_road_P -2.579e-05 9.662e-06 -2.669 0.007608 **
Dist_road_S -8.984e-05 6.592e-05 -1.363 0.172909
Dist_agri -4.340e-05 5.029e-05 -0.863 0.388227
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2607.7 on 1948 degrees of freedom
AIC: 2627.7

Number of Fisher Scoring iterations: 4

Cómo se observa, las variables que influyen significativamente sobre la deforestación en el periodo 1975-1985 son la distancia a carreteras principales (Dist_road_P) y, de manera más marginal (p = 0.081), la orientación (Aspect). La primera variable influye negativamente sobre la probabilidad de deforestacion (Estimate = -2.579 e-05), es decir, que a medida que nos alejamos de las carreteras principales la probabilidad de deforestación disminuye. La segunda es más difícil de intepretar, puesto que la orientación va de 0º (Norte) a 360º (Norte). En este caso, sería más interpretable utilizar una variable como insolación, que está directamente relacionada con la orientación y la pendiente.

Puesto que el resto de variables no son significativas se puede reducir el modelo.

model.red1 <- glm(Change ~ Aspect + Dist_road_P, data = fnf75_85, family = binomial)
summary(model.red1)

Call: glm(formula = Change ~ Aspect + Dist_road_P, family = binomial, data = fnf75_85) Deviance Residuals: Min 1Q Median 3Q Max -1.5055 -1.3468 0.9385 1.0064 1.2628 Coefficients: Estimate Std. Error z value Pr(>|z|)

(Intercept) 7.514e-01 1.092e-01 6.881 5.92e-12 ***
Aspect -7.502e-04 4.593e-04 -1.633 0.10236
Dist_road_P -2.496e-05 8.126e-06 -3.072 0.00213 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2615.3 on 1955 degrees of freedom
AIC: 2621.3

Number of Fisher Scoring iterations: 4

Y se comparan los dos modelos con la función anova().

anova(model.full, model.red1, test = "Chi")

Analysis of Deviance Table

Model 1: Change ~ Elevation + Slope + Aspect + Dist_river + Dist_City.20T +
Dist_village + Dist_road_P + Dist_road_S + Dist_agri
Model 2: Change ~ Aspect + Dist_road_P
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 1948 2607.66
2 1955 2615.32 -7 -7.66 0.36

Lo que indica que el modelo más simple no es significativamente peor que el modelo completo (p = 0.36). Así que aceptamos este modelo por el momento. Sin embargo, cuando analizamos los resultados se observa que la variable orientación no es significativa (p = 0.102), por lo que probamos a reducir el modelo un poco más.

model.red2 <- glm(Change ~ Dist_road_P, data = fnf75_85, family = binomial)
summary(model.red2)

Call: glm(formula = Change ~ Dist_road_P, family = binomial, data = fnf75_85) Deviance Residuals: Min 1Q Median 3Q Max -1.4540 -1.3500 0.9432 1.0034 1.2042 Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.333e-01 8.129e-02 7.791 6.66e-15 ***
Dist_road_P -2.551e-05 8.115e-06 -3.144 0.00167 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2618.0 on 1956 degrees of freedom
AIC: 2622

Number of Fisher Scoring iterations: 4

Y volvemos a comparar los dos modelos:

anova(model.red1, model.red2, test = "Chi")

Analysis of Deviance Table

Model 1: Change ~ Aspect + Dist_road_P
Model 2: Change ~ Dist_road_P
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 1955 2615.32
2 1956 2617.99 -1 -2.67 0.10

Así que de nuevo nos encontramos que el modelo reducido (model.red2) no es significativamente peor que el modelo anterior (model.red1). Aceptamos este último como válido, por lo que concluimos que la probabilidad de deforestación en el periodo 1975-1985 estuvo influida fundamentalmente por la presencia de carreteras principales (¿se debió la deforestación a un uso maderero por parte de empresas forestales?).

Por último, se representan gráficamente los resultados.

xv<-seq(0,30000,10)
yv<-predict(model.red2, list(Dist_road_P = xv), type = "response")
plot(as.integer(fnf75_85$Change) ~ fnf75_85$Dist_road_P, type = "p", ylab = "Probability of forest loss", xlab = "Distance to main roads (m)")
lines(xv,yv)

¿Qué nos dice este gráfico? Básicamente, que la probabilidad de deforestación disminuye con la distancia a carretera principales, pero disminuye poco. Así, cuando uno se aleja 25 kms de las carreteras principales (una distancia bastante grande), la probabilidad de deforestación disminuye del 65% a apenas el 50%.

Cómo hay muchos puntos en el gráfico, no se aprecia realmente si hay más 0s que 1s en distancias próximas a las carreteras principales y a la inversa... sería bueno intentar ajustar funciones de densidad de los 0s y 1s en el mismo gráfico.

¿Qué ocurre si la variable respuesta es una proporción?

Puede ser que la variable respuesta no sean 0s y 1s sino que sea una proporción o porcentaje, ej. porcentaje de semillas que sobreviven en un experimento, sex ratios, proporción de píxeles deforestados en un determinado área, etc. En este caso, tambińe podemos llevar a cabo un análisis glm con distribución de errrores de tipo binomial. Las únicas diferencias con respecto a los pasos descritos anteriormente serían:

Para definir la variable respuesta creamos un objeto con dos columnas, utilizando la función cbind() para juntar dos vectores que contengan el número de éxitos y fracasos en cada uno de los casos.

n.fracasos <- denominador - n.exitos
var.respuesta <- cbind(n.exitos, n.fracasos)

Comprobar que no hay sobredispersión (la devianza residual no es mayor que los grados de libertad residuales) y corregirlo en caso necesario utilizando una familia de tipo cuasibinomial (family = quasibinomial).

Bibliografía recomendada:

- Crawley, M.J. 2007. The R Book. Wiley, London.

Comparación de la pendiente de dos modelos (normal y poisson) de riqueza de especies

2009-01-19T19:37:00.019+01:00

En un estudio que está llevando a cabo Fabio Suzart de Alburquerque (Departamento de Ecología, Universidad de Alcalá) se quiere comparar la pendiente de la recta de regresión de la riqueza de especies nativas y exóticas en el Reino Unido frente a una variable explicativa relacionada con la entrada de energía en el sistema (el primer eje de un análisis de componente principales), correlacionado fundamentalmente con la temperatura media anual, la temperatura media del mes de enero, y la evapotranspiración potencial). El principal problema con el que se enfrenta es que la distribución de la riqueza de nativas sigue una distribución normal, mientras que la distribución de especies exóticas sigue una distribución más de tipo Poisson o log-normal (ver Figura 1).

El siguiente código de R permite la lectura de los datos (accesibles con conexión a internet), el ajuste de dos modelos "normales" para la riqueza de nativas y exóticas, y la creación de: (1) histogramas de los residuos de los modelos de riqueza de especies nativas y exóticas; y (2) gráficos de correlación de estas variables con el primer eje del análisis de componentes principales (PCA).

data <- read.table(url("http://archivos-para-subir.googlegroups.com/web/datos+regresion.txt?gda=x8KpMUUAAACkOfwqHAVd4YqgfIB09GDRRDpDrIn9MKnJm3pWH7Uz3yEZ0eBGq9HKqr7NTMdj8uts7Hm76flVp_vYI1lbK8kaGu1iLHeqhw4ZZRj3RjJ_-A"), header = T, sep = "\t")

glm.native <- glm(data$NAT_RICH ~ data$Factor1.energy) glm.exotic <- glm(data$ALIEN_RICH ~ data$Factor1.energy)
par(mfcol= c(2, 2))
hist(residuals(glm.native), main = "Histogram of native SR residuals", xlab = "Number of species")
par(new = T)
plot(density(residuals(glm.native)), xlab = "", ylab = "", main = "", axes = F)
hist(residuals(glm.exotic), main = "Histogram of exotic SR residuals", xlab = "Number of species")
par(new = T)
plot(density(residuals(glm.exotic)), xlab = "", ylab = "", main = "", axes = F)
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
plot(data$ALIEN_RICH ~ data$Factor1.energy, main = "Exotic species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")

Figura 1. Histogramas de los residuos de los modelos de especies nativas (arriba izquierda) y exóticas (abajo izquierda) y gráficos de correlación de la riqueza de especies nativas (arriba derecha) y exóticas (abajo derecha) con el primer eje del análisis de componentes principales (PCA).

Una posible solución a este problema es utilizar modelos lineales generalizados (GLM) con una distribución de errores de tipo Poisson y una función de vínculo de tipo logarítmica para linealizar la relación entre la riqueza de especies nativas y la variable independiente. Esto impide comparar la pendiente de ambas rectas de regresión de manera directa pero nos permite, en cambio, establecer hipótesis sobre los procesos que generan la respuesta de las variables dependientes (riqueza de nativas frente a riqueza de exóticas).

glm.native <- glm(data$NAT_RICH ~ data$Factor1.energy)
glm.exotic <- glm(data$ALIEN_RICH ~ data$Factor1.energy, family = poisson)

El primer modelo es un modelo normal. Hay que considerar que los GLM no son más que una extensión de los modelos lineales, por lo que los modelos lineales que comúnmente se utilizan (ANOVA, regresión) no son más que un caso particular de GLM en dónde la distribución de errores es de tipo gaussiana (normal). Con la función summary() se obtiene un resumen de los modelos ajustados.

summary(glm.native)

Call:
glm(formula = data$NAT_RICH ~ data$Factor1.energy)

Deviance Residuals:
Min 1Q Median 3Q Max
-428.470 -45.604 5.102 47.498 218.190

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 421.489 1.515 278.12 <2e-16>
data$Factor1.energy 59.039 1.516 38.95 <2e-16>
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 5172.153)

Null deviance: 19483365 on 2251 degrees of freedom
Residual deviance: 11637345 on 2250 degrees of freedom
AIC: 25652

Number of Fisher Scoring iterations: 2

summary(glm.exotic)

Call:
glm(formula = data$ALIEN_RICH ~ data$Factor1.energy, family = poisson)

Deviance Residuals:
Min 1Q Median 3Q Max
-18.7702 -4.5256 -0.3740 3.1827 23.9153

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.900846 0.001945 2519.6 <2e-16 ***
data$Factor1.energy 0.543437 0.002224 244.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 158010 on 2251 degrees of freedom
Residual deviance: 87121 on 2250 degrees of freedom
AIC: 101806

Number of Fisher Scoring iterations: 5

Como se puede observar, las pendientes de la recta de regresión de ambos modelos no son comparables, entre otras cosas, porque en el modelo Poisson la variable respuesta (riqueza de exóticas) está transformada logarítmicamente. Esto ya nos da una idea sobre el comportamiento de ambos tipos de riqueza de especies. Así, mientras que la riqueza de especies nativas varía linealmente con el primer eje del PCA (energía), la riqueza de exóticas cambia exponencialmente, es decir, que para valores bajos de la variable explicativa, el incremento de la riqueza será relativamente pequeño, pero para valores altos de la variable explicativa, el incremento de la riqueza de especies será muy alto.

Esto se puede observar en el modelo ajustado a los puntos para cada uno de los dos tipos de riqueza.

par(mfrow=c(1,2))
yv.native <- predict(glm.native, type= "response")
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.native))
yv.exotic <- predict(glm.exotic, type= "response")
plot(data$ALIEN_RICH ~ data$Factor1.energy, main = "Exotic species richness (response) vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.exotic))

Y alternativamente podemos dibujar las predicciones del modelo en escala logarítmica. Esta relación es lineal, pero hay que tener en cuenta que la variable respuesta (riqueza de especies) está en escala logarítmica, por lo que la comparación de las pendientes de ambas rectas no es posible.

par(mfrow=c(1,2))
yv.native <- predict(glm.native, type= "response")
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.native))
yv.exotic2 <- predict(glm.exotic, type= "link")
plot(log(data$ALIEN_RICH) ~ data$Factor1.energy, main = "Exotic species richness (link) vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "log(Number of species)")
lines(sort(data$Factor1.energy), sort(yv.exotic2))

Este estudio forma parte del siguiente trabajo: Suzart de Alburquerque, F., Castro, P. Cayuela, L. & Rodríguez, M.A. Assessing the influence of environmental and human-related factors on native and exotic species richness in Great Britain. In prep.