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

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

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

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

La PRCF ¿qué es y cómo se calcula?

Recientemente me he puesto a analizar unos datos de series temporales de dinámica de poblaciones de procesionaria (estimadas a partir de mediciones del grado de infestación de rodales) en Mora de Rubielos, Teruel, en el marco de una colaboración a tres bandas entre el Laboratorio de Sanidad Vegetal de Mora de Rubielos (Rodólfo Sánchez, Gerárdo Sánchez), la Universidad de Granada (José Antonio Hódar, Regino Zamora) y la Universidad Rey Juan Carlos (en mi humilde persona).

Para poder discernir entre factores endógenos y exógenos que podrían estar causando las oscilaciones de las series, una de las primeras cosas que hay que hacer es tratar de determinar si hay dependencia de la densidad y con cuantos retardos. Es decir, si el hecho de que, por ejemplo, uno, dos o tres años atrás, la población estuviera en un máximo, va a influir en que ahora, en tiempo t, la población haya sufrido un descenso en el número de individuos como consecuencia del agotamiento de los recursos u otros factores que podrían estar operando a nivel intraespecífico. Pues bien, para determinar con cuantos retardos la propia variable tiene un efecto sobre la observación de esta misma variable en tiempo t, Alan Berryman (1999) sugiere el uso de lo que llama la PRCF (partial rate correlation function), que no traduzco porque no se si lo haría bien. Para aquellos que estéis un poco familiarizados con la jerga de series temporales, no os será desconocido la función de autocorrelación (ACF) y la función de autocorrelación parcial (PACF). Estas son herramientas exploratorias básicas para identificar la estructura de una serie temporal. Pues bien, en el caso de dinámica de poblaciones, Berryman sugiere que es más apropiado utilizar la PRCF... ¿y qué es esto de la PRCF? Se trata de correlacionar la tasa de cambio poblacional en tiempo t, con la densidad poblacional en tiempo t-1, t-2 y así sucesivamente. La tasa de cambio de la población es lo que Berryman denomina como variable Rt, y que en realidad es el log(Nt/Nt-1). Es decir, una medida de cómo la población cambia entre un periodo y el siguiente.

En mi periplo analítico, me encontré con que no tenía ni idea de cómo calcular la PRCF en R, que es ya el único software con el que trabajo. Busqué en los foros y nada. Lo más que encontré fue la página de un software (PAS) que diseñó Alan Berryman ya hace tiempo y que explicaba con fórmulas cómo calcular la PACF y, a partir de ahí, la PRCF. Pero honestamente, no entendí mucho. Así que di en mis indagaciones con un artículo de Mauricio Lima publicado en PLoS ONE donde decía explícitamente que ellos habían calculado la PRCF con un código escrito en R. Escribí a Mauricio y este rápidamente me dirigió a Sergio Estay, quien con mucha diligencia me pasó el código, que con permiso suyo y agradeciéndole de antemano su enorme ayuda, pongo aquí a disposición de todos aquellos interesados (he modificado un par de cosas para generalizarlo, pero en esencia, todo el código es suyo).

Código para calcular la función PRCF en R 

Al final, y viendo el código de Sergio (y con alguna aclaración suya) entendí que era esto de la PRCF... y es bastante simple. Para más de un retardo la PRCF es equivalente a la PACF. Y para un retardo es simplemente la correlación entre Rt (la tasa de cambio) y Nt-1. La (última) pregunta que yo me hice fue ¿por qué para dos o más retardos la PRCF es igual a la PACF? Y dándole vueltas mientras volvía a casa en coche, llegué a una conclusión, que expongo a continuación y que no se si será más o menos acertada (si alguien conoce una respuesta más acertada por favor que me corrija).

Básicamente lo que se está diciendo es que para dos retardos (y generalizable a otros retardos por encima de uno):

r(Rt ~ Nt-2) habiendo quitado el efecto de Nt-1 sobre Rt 
r(Nt ~ Nt-2) habiendo quitado el efecto de Nt-1 sobre Nt 

¿Tiene esto algún sentido? Lo tiene si pensamos que Rt se calcula como el log(Nt/Nt-1). Por tanto si quitamos el efecto de Nt-1 sobre Rt, lo que nos queda es básicamente Nt. Por tanto, con dos o más retardos, la PRCF será equivalente a la PACF.

Y ya para terminar, incluyo aquí los resultados de la PRCF para las 11 parcelas con las que estamos trabajando en Mora de Rubielos, y en donde se ve que la dinámica de la procesionaria está regulada por efectos de densidad con un único retardo.

 

Curso de series temporales aplicadas a datos ambientales

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

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

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

Esperamos que este curso sea de vuestro interés.

Cómo estimar parámetros en un modelo por máxima verosimilitud (versión fuerza bruta)

Una parte importante de los modelos estadísticos es la forma en la que se estiman los parámetros. En modelos lineales (regresión, ANOVA, ANCOVA), la estimación se hace frecuentemente por el método de mínimos cuadrados. En modelos lineales generalizados (GLM) la estimación se hace por el método de máxima verosimilitud. Cómo opera este método es algo normalmente desconocido al usuario estándar de programas estadísticos (incluido R). Aunque hay muchas formas de estimar los parámetros de un modelo utilizando máxima verosimilitud, la idea subyacente es la misma: encontrar los parámetros que maximizan la probabilidad de los datos observados. La forma más intuitiva de hacer esta estimación es por fuerza bruta, esto es, probando muchas posibles combinaciones de los valores de los parámetros y ver qué combinación maximiza la probabilidad de los datos. A continuación pongo un breve manual sóbre cómo hacer esto en R con objeto de poder entender un poco mejor este proceso. 

Cómo estimar los parámetros de un modelo por máxima verosimilitud en R

No obstante, esta forma iterativa es computacionalmente inpracticable, como se explica en el manual, y distintos métodos se han ideado para solucionar estos problemas y encontrar parámetros óptimos (locales y globales). No entraré a detallar los distintos métodos de estimación por máxima verosimilitud, pero si entendemos la idea de fondo, creo que será mucho más fácil entender el proceso por el que generamos nuestros modelos y algunos de los problemas (por ejemplo, convergencia) que a veces surgen por el camino.

Modelos a la carta: El modelo gaussiano

Los modelos lineales generalizados (GLM) tienen tres componentes: (1) la estructura de errores; (2) el estimador lineal y; (3) la función de vínculo. La función de vínculo linealiza la relación entre la respuesta y el estimador lineal. La función de vínculo, junto con el estimador lineal, constituyen lo que lo algunos autores denominan el "modelo científico". La máxima versatilidad en la modelización se alcanza cuando uno ajusta modelos a la carta, fuera de las restricciones que los GLM imponen. Ello implica que debemos escribir nuestras propias funciones, para después ajustarlas a los datos. El proceso se completa en el contexto de la selección de modelos por criterios de información cuando comparamos toda una batería de modelos plausibles ajustados por métodos de máxima verosimilitud. En R esto puede hacerse por medio del paquete 'likelihood' escrito por Lora Murphy. El paquete todavía no está disponible en CRAN, aunque es de esperar que lo esté pronto. Mientras tanto, este paquete se puede descargar de la página del curso de "Modelos y métodos de máxima verosimilitud" de Charlie Canham. Este curso fue por cierto impartido en abril en Granada y he de decir que ha sido uno de los mejores cursos a los que he asistido... ¡gracias Charlie!

Bueno, a lo que iba... en lo que respecta a los modelos científicos. Podemos usar muchos tipos de modelos dependiendo de los datos que tengamos. No tienen por qué ser lineales. Una familia de modelos muy versátiles que permiten ajustar respuestas que aumentan primero y disminuyen después son los modelos gaussianos. La función gaussiana se expresa así:

$y=a\cdot e^{(\frac{-{(X-b)}^2}{(2\cdot c^2)})}$

Dicha función consta de 3 parámetros, que llamaremos 'a', 'b' y 'c', pero que podríamos haber denominado de cualquier otra forma. El parámetro 'a' indica el valor máximo que alcanza la curva. El parámetro 'b' determina la forma de la curva: monotónica, creciente, decreciente. El parámetro 'c' va a condicionar la pendiente de la curva (lo que en inglés se denomina sharpness).

Veamos ahora cómo cambios en los valores de estos parámetros determinan cambios en la forma de las curvas tipo. Tomemos un ejemplo en dónde vamos a modelar la abundancia de arácnidos en sistemas agrícolas (y) en función del tiempo (X). La idea con el siguiente código en R es que nos familiaricemos con esta función y sus parámetros.

Si modificamos el parámetro 'a' podemos obtener algunas las siguientes curvas:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=125, y=1025, "a=1000")
a <- 800
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=125, y=825, "a=800")
a <- 600
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=3)
text(x=125, y=625, "a=600")
a <- 400
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=125, y=425, "a=400")
a <- 200
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=125, y=225, "a=200")

Modificando el parámetro 'b' obtenemos cambios en la forma de la curva gaussiana como se observa en la siguiente gráfica:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=125, y=1025, "b=125")
b <- 200
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=200, y=1025, "b=200")
b <- 275
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=240, y=825, "b=275")
b <- 50
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=50, y=1025, "b=50")
b <- -25
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=15, y=825, "b=-25")

Cambios en el parámetro 'c' van a provocar cambios en las pendientes de la curva, como se observa a continuación:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=45, y=825, "c=0.01")
c <- 0.001
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=125, y=1025, "c=0.001")
c <- 0.25
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=125, y=150, "c=.025")
c <- 0.1
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=90, y=250, "c=0.1")
c <- 0.02
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=40, y=400, "c=0.02")

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

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:

1. 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.
2. 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).
3. Se calculan las distancias sobre este nuevo espacio geométrico y se calcula una matriz de distancia Y .
4. Se comparan las matrices de distancia X e Y y se mide cómo son de parecidas entre ellas (stress).
5. A partir de la con guración inicial, se reasignan los sitios (unidades muestrales) para reducir las distancias con la matriz X.
6. 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í.

Análisis de datos ecológicos en R

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

1. Una introducción a R (2.0 Mb).
2. Modelos lineales: Regresión, ANOVA y ANCOVA (621 Kb).
3. Modelos lineales generalizados (GLM) (416 Kb).
4. Introducción al diseño de experimentos (13.7 Mb).
5. Modelos lineales mixtos en R (1.9 Mb).
6. 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.

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

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

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