¿Me paso a la estadística Bayesiana?

A menudo encontramos trabajos en ecología que “aplican un marco bayesiano”, “corren cadenas de montecarlo”, o “infieren la heredabilidad con winbugs”. Cuando uno piensa si quiere o debe aplicar estadística bayesiana para algún problema en concreto y hace una búsqueda general, se encuentra con una auténtica batalla filosófica y metodológica entre unos bayesianos y otros frecuentistas. Los argumentos incluyen ideas profundas, a menudo discusiones sobre quién de los dos refleja realmente cómo funciona el método científico. Pero a menudo nos encontramos dos tipos de discusiones. Unas demasiado técnicas, inabordables para el público general, y otras demasiado generales que no nos permiten entender el problema ni dar respuesta a la pregunta de arriba. En esta entrada intento hacer algo intermedio: explicar los fundamentos para que los entienda alguien que habitualmente utiliza estadística (tradicional o frecuentista) y le pica la curiosidad cada vez que lee un paper que aplica el “método bayesiano”.

 

fig

https://xkcd.com/1132/

 Maximum likelihood vs distribuciones a posteriori

 En esencia, los entornos frecuentista y Bayesiano difieren en la forma de hacer inferencia, es decir, en cómo se calcula un parámetro poblacional como la media, la varianza, o la pendiente entre dos variables, a partir de una muestra. Imaginemos un modelo lineal simple, que viene dado por la ecuación de la recta y = a + bx + e (donde a y b son intercepto y pendiente, respectivamente, y e ~N(0, σ) es el error o residuo). Este modelo va a ser así utilicemos un entorno u otro, la diferencia será la forma en que inferimos los parámetros a, b y e a partir de los datos.

El entorno frecuentista establece que existe una única combinación de parámetros que minimiza las distancias de los datos a la recta, o lo que es lo mismo, que maximiza la verosimilitud (o likelihood). Es decir, siempre que intentemos trazar una línea a través de una nube de puntos, habrá una única línea que minimice las distancias a cada punto. Esta línea se define por un a, una b y dejará un e que es la distribución de distancias de los puntos a la línea, el residuo del modelo. Si el modelo cumple el requisito de tener residuos normales y homocedásticos (es decir, e es siempre una distribución normal de media 0 y desviación constante, σ), podremos generar algo muy parecido a nuestros datos a partir del modelo. En otras palabras, el modelo explica nuestros datos, es verosímil.

Lo anterior es simplemente la definición de verosimilitud (likelihood), utilizando la notación habitual: queremos explicar la variable y a partir de la variable x con un modelo lineal de parámetros a, b y e, que se suelen agrupar en el conjunto θ = {a, b, e}. Con ello, el likelihood se define como la probabilidad conjunta de los datos condicionada al modelo:

 p(y1, y2, …, yn |θ) = p(y1|θ) . p(y2|θ) … p(yn|θ)

Esto es, la probabilidad de que el dato yi pertenezca al modelo θ (multiplicamos todas esas probabilidades para obtener la probabilidad conjunta) Si vamos dando valores a los parámetros θ, veremos que para una determinada combinación de valores de a, b y e, la función alcanza un máximo. Si estamos en el entorno frecuentista, diremos que esos valores de a, b y e maximizan la verosimilitud. Habremos aplicado el método de Maximum Likelihood para ajustar el modelo.

Esta es la clave, estamos hablando de “la probabilidad de que el dato yi pertenezca al modelo θ”, p(y|θ), no estamos midiendo “la probabilidad de que el modelo se ajuste a los datos empíricos” p(θ|y). Pero, por qué íbamos a querer determinar la probabilidad de los datos, que hemos obtenido de forma empírica, a partir de un modelo que desconocemos? Por qué no hacer lo contrario, medir la probabilidad del modelo a partir de los datos que hemos obtenido en el campo?

Pensemos qué ocurre si lo que queremos es determinar en qué medida el modelo está soportado por los datos empíricos. Es decir, vamos a dar la vuelta a la tortilla del likelihood y medir la probabilidad del modelo condicionado a los datos, p(θ|y), en lugar de la probabilidad de los datos condicionados al modelo, p(y|θ). Cuando hablamos de probabilidades condicionales p(θ|y) ≠ p(y|θ), es decir, la probabilidad de que llueva un día nublado, p(lluvia | nubes), no es igual que la probabilidad de que esté nublado un día lluvioso p(nubes | lluvia). Del mismo modo, la probabilidad de los datos condicionados al modelo, p(y|θ) o likelihood, no es igual a la probabilidad del modelo condicionado a los datos, p(θ|y), que denominamos función a posteriori. La relación entre ambas probabilidades viene dada por el teorema de Bayes:

Donde p(θ) es la probabilidad del modelo a priori y p(y) es la probabilidad de los datos o evidencia. La probabilidad a priori mide la certidumbre previa acerca del modelo, la confianza que tenemos acerca de una cuestión general previamente a tener los datos empíricos. La evidencia es la probabilidad de los datos determinada a partir de la suma de todos los posibles valores de los parámetros multiplicados por la confianza previa acerca esos valores.

Sintetizando: en el entorno Bayesiano se está midiendo de forma explícita (1) si existe información previa acerca del tema en cuestión que se pueda comparar numéricamente con la información nueva que se aporta, p(θ); (2) cuál es la verosimilitud de nuestros datos en el contexto de esa información previa, p(y|θ); y (3) si nuestros datos tienen fuerza para falsar el conocimiento previo sobre el tema, p(θ|y).

Pero entonces, me paso a la estadística Bayesiana?

Casi siempre que surge la pregunta de si es mejor aplicar un método u otro para analizar unos datos, la respuesta es la misma: si los resultados no son los mismos aplicando cualquiera de los métodos, sospecha de los datos, no de los métodos. Seguramente el problema es que los datos no contienen demasiada información, no hay un tamaño muestral adecuado, o incluso la pregunta no es la adecuada. Si la pregunta es la adecuada y los datos tienen la respuesta, el método Bayesiano y el frecuentista deberían darte un resultado muy parecido, aunque cada uno tenga su manera de hacerlo.

Dicho esto, realmente podemos sacar bastante provecho del método Bayesiano para abordar preguntas que no son tan fáciles de resolver mediante estadística tradicional o frecuentista. El mayor inconveniente que tiene el entorno Bayesiano, es que está mucho menos trillado que la estadística tradicional, los paquetes de software requieren información sobre cosas que nos cuesta definir (como cuál es la distribución adecuada para mi prior, o cómo parametrizo la distribución gamma), y las ayudas de los paquetes estadísticos no son tan manejables e intuitivas como las de los paquetes tradicionales. Aun así, a veces interesa bucear en cuestiones sobre estadística Bayesiana, sobre todo si te interesa investigar la incertidumbre de los modelos, o los componentes de la varianza en una muestra. Esta es la receta que puedo dar, según mi experiencia:

 

  1. Si las técnicas que utilizas habitualmente son modelos lineales o modelos mixtos donde sólo te interesan los efectos fijos, no pierdas mucho tiempo en aprender estadística Bayesiana. En este caso, ajustar tus modelos con estadística Bayesiana no te va a aportar respuestas que no te aporte ya la estadística tradicional.
  2. Si en tu trabajo exploras en detalle los factores aleatorios de los modelos mixtos, e.g. analizas la estructura de la varianza usando modelos jerárquicos o modelos animales con pedigrees, dedica un tiempo a pensar sobre estadística Bayesiana, porque te puede aportar ideas interesantes. Por ejemplo, paquetes de R como MCMCglmm (Hadfield 2010) calculan modelos mixtos aplicando métodos Bayesianos y te dan toda la libertad para ajustar la estructura de la matriz de varianza-covarianza. Esto significa que puedes testar cuánto explican tus factores aleatorios y cómo son sus interacciones, tanto entre factores aleatorios como de los aleatorios con los fijos.
  3. Si utilizas modelos no lineales para analizar crecimiento de estructuras biológicas o poblaciones, por ejemplo, funciones logísticas, funciones de Bernoulli o Gompertz, entra en profundidad en estadística Bayesiana para encontrar todo un mundo en la selección de modelos (con programas como JAGS o WinBUGS). Por ejemplo, aplicando estadística Bayesiana y selección de modelos con criterios de información como el DIC, puedes testar la complejidad de los modelos no lineales, interacciones entre parámetros, y validar el ajuste de unos u otros modelos de una manera más sofisticada que con métodos frecuentistas tradicionales, que te puede abrir puertas a nuevas preguntas.
  4. Por último, aprende en profundidad estadística Bayesiana si te interesa la modelización de procesos estocásticos. Por ejemplo, para analizar modelos matriciales de poblaciones o modelos de invasión de mutantes en dinámicas evolutivas donde los parámetros son estocásticos. A menudo estos modelos se parametrizan con datos reales, por ejemplo, valores reales de tasas de natalidad, mortalidad, migración o mutación de genes. Si para una tasa demográfica puedes obtener, no un dato, sino una distribución de datos posibles, puedes parametrizar el modelo, no con un dato, sino con toda una distribución. Con eso vas a tener un output que no sólo describe la trayectoria de la población, sino su varianza en base a la dispersión de los datos empíricos con los que has parametrizado el modelo.

 

 

 

 

Anuncios

Cosas que conviene saber al usar AIC, DIC y otros criterios de información

Los indicadores con criterio de información abordan el compromiso entre la complejidad y la capacidad predictiva de un modelo. Cuanto más complejo sea un modelo peor será su capacidad para predecir en un abanico amplio de situaciones. Es decir, cuantas más variables y cuantas más interrelaciones entre los componentes incorpore un modelo más concretas serán sus predicciones pero también menos generalizables. Si un modelo es sencillo, e incluye sólo los componentes con mayor importancia en el sistema, podrá predecir en un mayor número de escenarios aceptando un cierto error de precisión.

Este compromiso es la base de que busquemos simplificar al máximo los modelos estadísticos reduciendo el número de variables predictoras. Eliminaremos variables del modelo con un ojo siempre puesto en su capacidad para explicar los datos. Cuando eliminando variables empecemos a perder capacidad explicativa, entenderemos que hemos alcanzado el mejor modelo posible, el más parsimonioso, aquel que explica lo máximo con el menor número de parámetros.

Para encontrar el modelo más parsimonioso de entre una batería de ellos necesitamos un indicador capaz de (1) medir la capacidad explicativa de un modelo y (2)  penalizar por su grado de complejidad. Esto es lo que hacen los indicadores con criterio de información (como AIC, BIC o DIC entre otros), que van a tener siempre la forma:

xIC = complejidad – bondad de ajuste

El objetivo es encontrar el modelo que hace mejores predicciones. Puesto que generalmente no se dispone de datos independientes para validar el modelo, podemos aproximar la calidad de las predicciones sobre los propios datos que se han utilizado para construirlo. Por lo tanto necesitamos penalizar el doble uso de los datos incorporando la medida de la complejidad. En general, los indicadores que existen difieren en la forma de medir la bondad de ajuste y la forma de penalizar por la complejidad.

Vamos a comparar AIC (Akaike Information Criterion), BIC (Bayesian Information Criterion) y DIC (Deviance Information Criterion). Los dos primeros utilizan Maximum Likelihood como criterio de bondad de ajuste y el número de parámetros como medida de complejidad. Es decir, se toman valores puntuales de los parámetros y de la bondad de ajuste que son los que maximizan el likelihood. El DIC, sin embargo, trabaja con la función completa de los parámetros, que construye mediante cadenas de Monte Carlo (MCMC). Esto permite medir la complejidad como el número y grado de interrelación entre los parámetros. El DIC es una generalización de los dos anteriores e incorpora una medida mucho más completa del ajuste y la complejidad.

Criterio de Información de Akaike

Fue el primero y sigue siendo el más utilizado en análisis de datos ecológicos. Viene dado por:

AIC = 2 k – 2 ln(L)

Donde k es el número de parámetros y L el valor máximo en la función de verosimilitud (Maximum likelihood). Por lo tanto AIC mide la bondad de ajuste a partir de la máxima verosimilitud del modelo, y la complejidad a partir el número de parámetros. La verosimilitud, p(y|θk), es la probabilidad de los datos condicionada a los k parámetros del modelo. Esta es una medida inversa a la suma de las distancias euclídeas de cada uno de los datos al modelo. Existe un único conjunto de valores de los parámetros que maximiza la función de verosimilitud (= que minimiza la suma de distancias euclídeas de los datos al modelo). Por lo tanto, cuanto más alto sea el valor de máxima verosimilitud mejor ajustado estará el modelo a los datos. Esto depende tanto de la forma del modelo como del grado de dispersión de los datos.

La medida de bondad de ajuste en este caso es 2 ln(L). Esta es una función creciente (la bondad de ajuste aumenta de forma logarítmica respecto al valor de máxima verosimilitud, L), pero la hacemos decreciente al multiplicarle -1. Así, cuando L tiende a infinito, la expresión -2 ln(L) que aparece en la función de AIC tiende a un valor muy pequeño (el indicador AIC se reduce cuando aumenta la bondad de ajuste). El valor de máxima verosimilitud, L, es el producto de las probabilidades de cada dato condicionado al modelo, por lo que L se obtiene multiplicando n valores entre 0 y 1 (o sumando sus logaritmos). Por lo tanto, destaco aquí para retomarlo después, que L depende tanto de la distancia de cada dato al modelo como del número de datos n.

La complejidad en el AIC viene dada por k, que es el número de parámetros del modelo. La penalización de 2k del AIC es equivalente a hacer la validación cruzada del modelo dejando un dato fuera (leave one out-cross validation). ¿Qué ocurre si aumentamos el número de parámetros, k, y eso no se traduce en un incremento de la verosimilitud, L? Es decir, si hacemos más complicado el modelo sin que con ello consigamos que explique mejor nuestros datos. Pues que el AIC aumentará del orden de 2 ∆k. Por tanto, diferencias de AIC de aproximadamente 2 entre modelos que difieren en un solo parámetro indican que la verosimilitud del modelo no cambia a pesar de que añadamos o quitemos un parámetro. Si lo que estamos haciendo es quitar parámetros para simplificar el modelo, podemos establecer como criterio dejar de quitar variables cuando el decremento de AIC < 2, lo que indica cambios no significativos en la bondad de ajuste del modelo.

Como veremos, AIC tiene una forma menos restrictiva que otros indicadores de medir la complejidad del modelo: sumando al indicador un número entero k que corresponde al número de parámetros. El problema de utilizar como única medida de complejidad es que su efecto es minúsculo comparado con el valor que toma -2ln(L) que depende del número de datos n (como he dicho antes L es proporcional a n). En otras palabras, perdonamos que el modelo sea muy complejo siempre y cuando tenga muchos datos para avalar cada parámetro. Esto se debe a que el AIC tiene como objetivo seleccionar el modelo que mejor hace predicciones dentro de los datos que lo han generado. Es decir, deja en el aire la cuestión de que un modelo complejo es también muy poco general, a pesar de que haya muchos datos avalándolo, y se centra en encontrar el modelo que hace mejores predicciones a pequeña escala. Tenemos que conseguir, por lo tanto, que aunque L dependa de n, esto no compense el exceso de complejidad del modelo.

Criterio de Información Bayesiano

Muy similar al AIC es el criterio bayesiano de información, BIC (Bayesian Information Criterion) de Schwarz, que viene dado por:

BIC = k ln(n) – 2 ln(L)

Donde, de nuevo k es el número de parámetros, L es el valor de máxima verosimilitud y n es el número de datos. Igual que el AIC se basa en la máxima verosimilitud como método de medida de la bondad de ajuste. Aquí vemos que la medida de la complejidad incorpora tanto k como ln(n). Esto independiza al indicador del tamaño muestral y hace que penalice más la complejidad que el AIC.

El resultado es que el BIC generalmente selecciona el modelo más abstracto, más sencillo y que hace predicciones a menor detalle, mientras que el AIC dará con un modelo más complejo y pragmático que hace predicciones con mayor detalle pero dentro de nuestros propios datos.

A pesar de ello, tanto AIC como BIC incurren en una excesiva simplificación del significado de complejidad del modelo. Recordemos que ambos miden la complejidad a partir del número de parámetros, k. Uno de los problemas de utilizar ese criterio de penalización es que require asumir que los parámetros son independientes entre ellos. Esto sólo es cierto si la colinealidad entre las variables del modelo es cero y además cada variable contribuye del mismo modo al modelo, lo cual no es cierto si se trata de modelos jerárquicos, donde la varianza está compartimentada por factores aleatorios.

Otra cuestión importante de medir la complejidad con el número de parámetros es que se ignora la incertidumbre en el modelo. Es decir, con AIC y BIC, el modelo se ha generado eligiendo el conjunto de valores de los parámetros que maximizan la verosimilitud, esto es, existe un valor puntual para cada parámetro. Los parámetros pueden no entenderse como valores puntuales sino como funciones de distribución que dependen de la dispersión de los datos o incertidumbre del modelo.

Criterio de Información de la Devianza

El DIC (Deviance Information Criterion) de Spiegelhalter es una versión generalizada del AIC y del BIC que incorpora la batería de procedimientos y la lógica de la estadística bayesiana. El DIC utiliza cadenas de Monte Carlo para buscar la distribución de los parámetros. Es decir, no se buscan los valores puntuales de los parámetros que maximizan la función de verosimilitud, sino que se trabaja con la función completa de verosimilitud. Así, recordamos que la verosimilitud es la probabilidad condicionada de los datos al modelo, p(y|θk), y que AIC y BIC miden la bondad de ajuste como la altura del punto de máxima verosimilitud, L. Por su parte, el DIC extrae la devianza, -2 ln[p(y|θk)], la cual tiene en cuenta la función completa de versomilitud. La medida de bondad de ajuste es el promedio de la devianza:

devianza

Para medir la complejidad del modelo, el DIC utiliza de nuevo la devianza de las funciones de los parámetros. Así, cuando los parámetros de un modelo están correlacionados, la búsqueda con MCMC genera distribuciones con alta devianza. Esto es, si los parámetros covarían habrá mayor variabilidad en el espacio de parámetros. Para penalizar la complejidad, el DIC calcula el número efectivo de parámetros, pD, que se define como el promedio de la devianza menos la devianza evaluada en el promedio de los parámetros:

pD

Con ello, el DIC se define como:

DIC

El número efectivo de parámetros tiene en cuenta simultáneamente el tamaño muestral, la relación (covariación) entre los parámetros, el número de parámetros y el tamaño de los efectos de las variables.

Por supuesto el DIC no está exento de críticas. En concreto, el sistema por el que penaliza el número de parámetros no equivale a una forma de cross-validation como ocurre en el AIC, lo que lleva a seleccionar modelos sobre-parametrizados. La forma de corregir esto es multiplicar por un factor el número efectivo de parámetros, tal que el indicador converja al método de leave one out al igual que el AIC. Existen múltiples revisiones del método implementadas en winBUGS y JAGS.

En resumen, el DIC se puede entender como una generalización de AIC y BIC que relaja la asunciones de independencia y certidumbre de los parámetros. El DIC hace una búsqueda del modelo verdadero implementando técnicas de MCMC que aportan una visión más completa de las propiedades de los modelos estadísticos.

Explorando el caos en la ecuación logística con R

Para abrir boca entramos en la resolución numérica del modelo de Verhulst, que describe el crecimiento logístico de una población cuya tasa de crecimiento disminuye al acercarse a un límite de capacidad de carga del medio. Aunque la función logística se puede resolver por vía analítica, el análisis mediante iteraciones nos deja explorar y jugar con la dinámica del modelo. Aunque me baso en las notas de Geoffrey R. Goodson (http://pages.towson.edu/goodson/Chaos_Notes.pdf) utilizo un código de R propio, que se puede descargar de github.

El mapa logístico es una familia de sistemas dinámicos no lineales de la forma µx (1 – x). El modelo de Verhulst, que viene dado por ry (1 – y / K), representa una población cuya dinámica está determinada por dos procesos: aceleración y ralentización de la tasa de crecimiento. Así, y(t) aumenta proporcionalmente al número de efectivos en cada instante y se ralentiza conforme se aproxima a la capacidad de carga del medio, K. Por lo tanto, podemos saber el número de individuos en un instante t + 1 a partir del instante anterior t, una tasa de crecimiento r y una K, que fijaremos en 1, trabajando por tanto con la ecuación generalizada.

Por ejemplo, una población que crece con una r = 1.5, partiendo de unas condiciones iniciales y(0) = 0.2; en t = 1 la población habrá aumentado a: y(1) = ry(0) (1 – y (0)/ K), es decir: y(1) = 1.5 * 0.2 * (1 – 0.2 / 1) = 0.24. Si seguimos iterando, en t = 2, la población llegará a y(2) =  ry(1) (1 – y (1)/ K), es decir, 1.5 * 0.24 * (1 – 0.24) = 0.27. Podemos seguir, y llegaremos a y(3) = 0.30; y(4) =  0.29; y(5) = 0.29 … y(inf) = 0.29. Es decir, la población crece desde y(0) = 0.2 a ritmo de r = 1.5, tras un pequeño reajuste en t = 3, rápidamente se estabiliza en 0.29.

Ese parece el típico caso de crecimiento logístico de una población. Veamos que ocurre para una tasa de crecimiento mucho mayor, r = 3.9. Igualmente partimos de y(0) = 0.2, por lo que y(1) = 3.9 * 0.2 * (1 – 0.2) = 0.62. La población ha llegado mucho más lejos en un sólo paso, lógicamente. Entonces y(2) = 3.9 * 0.62 * (1 – 0.62) = 0.92. Muy cerca del valor que le hemos dado a K, la población debería estabilziarse a partir de aquí. A ver, y(3) = 3.9 * 0.92 * (1 – 0.92) = 0.29. Pues no, hemos caído más abajo, y(4) = 3.9 * 0.29 * (1 – 0.29) = 0.80. y(5) = 0.62; y(6) = 0.92… y(inf) = ?? La población ha tomado una dinámica caótica. ¿Podemos predecir el futuro de la población? Esta es una situación de caos determinista.

Habitualmente se analiza el comportamiento del modelo logístico normalizado x ∈ [0,1] y para K = 1. Para valores de r ∈ (0, 4] es posible ver el espectro de dinámicas del sistema que van desde un único punto fijo que actúa como atractor (0 < r < 1), dos atractores fijos (1 < r < 2) y una dinámica cíclica (2 < r < 3) que se convierte en caótica (3 < r < 4). Vamos resolver la derivada de la ecuación logística de Verhulst a partir de la ecuación iterativa: y(t) = ry(t – 1) (1 – y(t – 1) / K). Esto significa que calcularemos y(1) a partir del valor que la función toma en el instante anterior y(0), que constituye la condición inicial, y continuar iterando y(2) = f(y(1)), y(3) = f(y(2)), etc. Pero antes de iterar directamente la función analizaremos su dinámica a partir de su representación geométrica.

Podemos representar la función y(t), que corresponde a la derivada de la función sigmoidea de crecimiento poblacional, y es por tanto una parábola que en nuestro caso pasa por 0 cuando y = 0 y cuando y = 1 y tiene un máximo en y = 1/2. Así, a partir de un valor inicial y(0), que encontramos en el eje de abscisas, podemos leer su correspondiente ordenada y(1) = f(y(0)) directamente a partir de la parábola. Del mismo modo, ahora podemos proyectar el valor de y(1) en la parábola opuesta (que se apoya en el eje de ordenadas) y buscar y(2) = f(y(1)) en el eje de abscisas. Si repetimos el proceso, nos nos conducirá al atractor:

# Empezamos con una r = 2, con la típica forma logística
r = 2 
K = 1
# Definimos y(t) como una secuencia en el intervalo [0, 1]
yt <- seq(0,1,0.001) 
# Y calculamos y(t + 1) a partir de la ecuación iterativa
yt1 <- r* yt *(1-( yt /K))
plot(yt1 ~ yt, ylim = c(0,1), xlim = c(0,1), ylab = "y(t + 1)", xlab = "y(t)", type = "l")

lines(yt ~ yt1, col = "blue")

Asignamos un valor inicial y0 y proyectamos para buscar el correspondiente y1:

y0 = 0.2
y1 = yt1[yt == y0]
arrows(y0, 0, y0, y1, col = "red")

Si damos la vuelta a la gráfica podemos buscar y2 a partir de y1:

# Busca el valor en la parábola más cercano a y1 
nexty <- which.min(abs(yt - y1)) 
y2 = yt1[nexty]
arrows(y0, y1, y2, y1, col = "red")

 

Fig1

 

Repitiendo el proceso se puede ver que llegaremos al punto en el que ambas parábolas intersectan. Esto es, los puntos de intersección entre las parábolas representan puntos de atracción, que en este caso (para r = 2) se encuentran en c = 0 y c = 1 – 1 / r = 0.5. Pero hay una forma más sencilla de cambiar entre ambas parábolas, que es representando la línea diagonal tal que y(t)= y(t + 1). Igual que antes, los puntos de intersección entre la parábola y(t) y la recta y(t) = y(t + 1) se corresponden con los atractores del sistema. Así:

plot(yt1 ~ yt, ylim = c(0,1), xlim = c(0,1), ylab = "y(t + 1)", xlab = "y(t)", type = "l")
l = yt
lines(l ~ yt)

Buscamos de nuevo y1 = f(y0):

y0 = 0.2
y1 = yt1[yt == y0]
segments(y0, 0, y0, y1, col = "red")
segments(y0, y1, y1, y1, col = "red")
prevy = y1

Y lo replicamos en bucle:

iteraciones = 50
for(i in 1:iteraciones){ 
 y1 <- which.min(abs(yt - prevy)) 
 nexty = yt1[y1] 
 segments(prevy, prevy, prevy, nexty, col = "red") 
 segments(prevy, nexty, nexty, nexty, col = "red") prevy = nexty
}
Rplot

Vemos que cuando 0 < r <1 existe un único atractor dentro de [0, 1] en c = 0. Es decir, cualquier punto de partida y(0) llevará a la población al colapso cuando la tasa de crecimiento es < 1. Sin embargo, cuando r aumenta por encima de 1 se crea un nuevo punto de atracción (las parábolas se cortan en dos puntos: c = 0; c = 1 – 1/r) en el que la población tenderá a estabilizarse. Por lo tanto, la población llegará a ese punto creciendo si y(0) < c o decreciendo si y(0) > c.

La cosa se complica cuando r es mayor que 2, donde aparece un punto de atracción cíclico. Esto es, el sistema se aproxima lentamente al punto de corte y, o bien converge al mismo (como ocurre para r = 2.9), o bien permanece fluctuando alrededor. Así, para 2 < r < 2.5 aproximadamente, el punto de corte atrae y repele simultáneamente al sistema haciéndolo orbitar alrededor del mismo. La población alcanzaría el punto de atracción y empezaría a producir ondas cada vez mayores hasta alcanzar un equilibrio fluctuante alrededor el punto de atracción. ¿Qué ocurre cuando r > 3.5? El efecto de repulsión del punto de corte aumenta generando cada vez ciclos más amplios y complejos, hasta que r cercanos a 4 generan una dinámica caótica en la población.

Vamos a iterar la función.

# Seleccionamos el valor de r y las condiciones iniciales y0
r = 0.9
y0 = 0.1
K = 1
# Generamos el vector que almacenará las iteraciones y le ponemos el y0 en la primera posición
y = numeric(100)
y[1] = y0
# Calculamos las y(t + 1) = f(y(t)) en bucle
for(i in 2:100) {
 y[i] = r*y[i-1]*(1-(y[i-1]/K))
}

plot(y, ylim = c(0,1), ylab = "y", xlab = "Tiempo, t", type = "l")

Fig4

La situación para r = 3.9 representa el caos determinista en la población. Determinista porque la dinámica del sistema puede conocerse y predecirse, en teoría, siempre que se conozcan las condiciones iniciales. Sin embargo, un pequeño cambio en esas condiciones cambiará todo el sistema en forma de “efecto mariposa”.

De forma muy sencilla podemos visualizar todos los puntos atractores del sistema simultáneamente para múltiples valores de r en el intervalo (1, 4]. Esto es, podemos representar el valor de y(t) que actúa como atractor para toda una serie de valores de r, de forma que 1 < r < 3 tendríamos un atractor en c = 1 – 1/r (sabemos que para 1 < r < 2 el sistema converge a un punto, a pesar de que pasa por un período de oscilaciones); cuando 3 < r < 3.5, aproximadamente, el sistema oscila y finalmente alcanza el caos 3.5 < r < 4.

# Generamos una función que itera la ecuación logística
logistica<-function(vector,y0,K,r){
 R=length(vector)
 rn=length(r)
 matriz=matrix(nrow=R,ncol=rn)
 matriz[1,]<-y0
 iterar<-function(a){
  r*a*(1-(a/K))}
  for(i in 1:R){
 if(i>1){ 
  matriz[i,]=iterar(matriz[i-1,])
 }else{
  matriz[i,]=iterar(matriz[i,])}
 }
 matriz
}
# Y otra función que aplica la anterior para cada valor de r y 
# representa el valor de y(t) que actua como atractor
bifurcacion=function(modelo){
 bifur=data.frame()
 L=length(modelo[,1])
 tabla=modelo[(L/5):L,]
 for(i in 1:length(r)){
 mat=data.frame(y=tabla[,i],r=r[i])
 bifur=rbind(bifur,mat)
 }
 plot(bifur,pch=".",cex=1)
}

 

# Generamos una secuencia de valores de r
r <- seq(1,4,0.005)
# Iteramos la logística con la función anterior
yt <- seq(0,1,0.001)
modelo=logistica(yt, y0 = 0.2, K = 1, r)
# Y generamos el diagrma de bifurcación
bifurcacion(modelo)


Fig5

 

 

 

 

 

 

 

 

 

¿Por qué un blog sobre modelos biológicos?

El carácter eminentemente empírico de la biología que nos enseñan en institutos y universidades ha alejado a los biólogos del enorme potencial teórico que tiene esta ciencia. Autores mucho más cercanos al razonamiento deductivo como Erwin Schrödinger (en Qué es la vida, 1944), Gregory Chaitin (en Proving Darwin 2012) o el propio Ernst Mayr (en This is biology 1997) se refieren a la “inmadurez” que a menudo parece mostrar una biología excesivamente centrada en detalles, mecanismos y las causas próximas. Pero la biología dispone de un enorme potencial para generar hipótesis por vía deductiva, y que se sustenta en la Teoría Evolutiva. Un buen ejemplo de ello es la obra de Richard Dawkins (como El Gen Egoísta 1976), donde queda patente la capacidad de la modelización matemática y la deducción para abordar las causas últimas (los por qués) de problemas como la evolución del altruismo, la evolución del sexo y el conflicto sexual o el conflicto entre hermanos.

A este blog se traerán soluciones matemáticas a problemas biológicos, teóricos y prácticos. Las entradas incorporarán funciones en R que en algún momento he preparado para abordar este tipo de cuestiones.