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

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

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?

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:

1. 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.
2. 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.
3. 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.
4. 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

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

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

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

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

Principales diferencias de Lyx con otros editores de texto

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

Cómo instalar Lyx

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

Lyx y Sweave

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

Instalar Sweave

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

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

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

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

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

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

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

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

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

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

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

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

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

data(Amazonia)
data(soils)

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

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

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

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

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

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

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

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

plot(Amazonia.rda)

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

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

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

Nuevo artículo publicado en Ecology and Society

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

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

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:

1. Análisis de componentes principales
2. Análisis de la varianza multivariado (MANOVA)
3. Escalamiento multidimensional no métrico (NMDS)
4. Análisis de correspondencias canónico (CCA)
5. Á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

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

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:

1. Conceptos estadísticos básicos: ANOVA y regresión
   
2. Cosas importantes antes de empezar
   
3. Cómo ajustar un modelo lineal en R
   
1. Un ejemplo de regresión
2. Un ejemplo de ANOVA
   
3. Un ejemplo de ANCOVA
   
4. Interacción entre factores o factores y co-variables
   
4. Evaluación de los supuestos del modelo: Exploración de los residuos
5. Problemas de colinealidad: Reducción de variables
6. Estandarización de coeficientes

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

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:

1. Conceptos básicos
2. La organización de R
1. Tipos de sistemas gráficos
2. Formato de gráficos
   
3. Gráficos tradicionales
1. El paquete graphics
2. Funciones de alto nivel: representación de una o dos variables
3. Funciones de alto nivel: representación de múltiples variables
4. Funciones de bajo nivel
5. Cómo personalizar un gráfico
   
4. Los paquetes grid y lattice

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

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:

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

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

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

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

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

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

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

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

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

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

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

Referencias:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Representamos ahora la matriz de correlaciones.


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

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

Como utilizar el potencial de GRASS desde R

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.