Por: Raymond L. Tremblay
Los modelos poblacionales bayesianos permiten integrar explícitamente la incertidumbre demográfica, de muestreo y de proceso en el análisis de la dinámica poblacional. Este capítulo introduce los modelos de proyección poblacional bayesianos (Bayesian MPP) como una extensión coherente de los modelos matriciales clásicos.
Los modelos de proyección poblacional bayesianos (Bayesian Population Projection Models, MPP) extienden los modelos matriciales tradicionales al tratar los parámetros vitales (supervivencia, crecimiento y fecundidad) como variables aleatorias en lugar de valores puntuales fijos. Esta formulación permite propagar explícitamente la incertidumbre desde los datos primarios hasta las proyecciones de abundancia y crecimiento poblacional (λ), algo fundamental en poblaciones con tamaños de muestra reducidos, series temporales cortas o alta variabilidad ambiental.
En este capítulo se presenta el marco conceptual de los Bayesian MPP, mostrando cómo se combinan modelos de observación y proceso dentro de una jerarquía bayesiana para estimar matrices poblacionales con distribuciones posteriores completas. Se discuten las ventajas clave de este enfoque frente a métodos frecuentistas, incluyendo la incorporación de información previa, el manejo coherente de datos faltantes y la cuantificación directa de la incertidumbre en predicciones a corto y largo plazo.
A través de ejemplos aplicados, se ilustra cómo los Bayesian MPP facilitan la evaluación del riesgo poblacional, la comparación entre escenarios de manejo y la integración posterior con análisis de elasticidad, sensibilidad y LTRE. Este capítulo se diferencia de los anteriores porque desplaza el énfasis desde la estimación puntual de parámetros hacia la inferencia probabilística completa de la dinámica poblacional, un enfoque especialmente relevante en biología de la conservación y ecología aplicada.
if (!require("pacman")) install.packages("pacman")
pacman::p_load(janitor, tidyverse)
library(tidyverse)
library(janitor)
library(popbio) # para la función projection.matrix()
library(raretrans) # disponible en CRAN: install.packages("raretrans")
library(DiagrammeR)
library(Rage)rlt_theme <- theme(
axis.title.y = element_text(colour = "grey20", size = 10, face = "bold"),
axis.text.x = element_text(colour = "grey20", size = 10, face = "bold"),
axis.text.y = element_text(colour = "grey20", size = 15, face = "bold"),
axis.title.x = element_text(colour = "grey20", size = 15, face = "bold")
) +
theme(
# Remover los bordes
panel.border = element_blank(),
# Remover las líneas
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
# Remover el fondo del panel
panel.background = element_blank(),
# añadir lineas más gruesas
axis.line.x = element_line(colour = "black", linewidth = 1),
axis.line.y = element_line(colour = "black", linewidth = 1)
)
diverging= c("#009392","#CF597E","#E9E29C", "#39B185","#EEB479" , "#9CCB86" )
# --- 2. Combinarlos en un objeto de tipo lista ---
# Esta lista se puede añadir al gráfico con un solo '+'
rlt_style_fill <- function() {
list(
rlt_theme,
scale_fill_manual(values = diverging)
)
}
rlt_style_colour <- function() {
list(
rlt_theme,
scale_colour_manual(values = diverging)
)
}En estadística clásica (frecuentista), tratamos los parámetros (e.g., la probabilidad de transición de plántula a juvenil) como valores fijos pero desconocidos, y los datos como variables aleatorias. La inferencia consiste en estimar el parámetro a partir de los datos: si observamos 2 transiciones de 7 plántulas, el estimador puntual es 2/7 = 0.286, sin más información.
En el marco bayesiano, invertimos la perspectiva: los datos son lo que observamos, y los parámetros se tratan como variables aleatorias sobre las que tenemos creencias. Esas creencias se expresan como una distribución previa (lo que sabíamos antes de mirar los datos) y se actualizan a una distribución posterior (lo que sabemos después). La regla de actualización es la fórmula de Bayes:
\[\underbrace{P(\theta \mid \text{datos})}_{\text{posterior}} \propto \underbrace{P(\text{datos} \mid \theta)}_{\text{verosimilitud}} \times \underbrace{P(\theta)}_{\text{previa}}\]
En español plano: la posterior es proporcional a la previa multiplicada por la verosimilitud. La previa aporta información cuando los datos son escasos; la verosimilitud (los datos) domina cuando hay muchas observaciones. Esta propiedad es exactamente lo que necesitamos para especies raras: si observamos 2 plántulas y 0 transiciones a juvenil, el estimador frecuentista da \(\hat{p} = 0\) (lo cual implica imposibilidad biológica), mientras que el estimador bayesiano combina las 2 observaciones con la previa y da una probabilidad pequeña pero positiva, con incertidumbre cuantificada.
La regla de Bayes en lenguaje cotidiano.
\[\text{Lo que crees ahora} \;\propto\; \text{Lo que creías antes} \;\times\; \text{Lo que dicen los datos}\]
Cuando los datos son muchos y consistentes, dominan: la posterior se parece a la verosimilitud. Cuando los datos son escasos o ambiguos, la previa importa: la posterior se parece a la previa. El método bayesiano hace este equilibrio explícito y reproducible.
Por qué leer este capítulo (incluso si nunca has hecho un análisis bayesiano). La estimación clásica de matrices de transición —projection.matrix()— es perfectamente adecuada cuando el tamaño de muestra es grande y cada transición se observó muchas veces. Pero en estudios reales de orquídeas, casi nunca es así: las muestras son pequeñas, varias transiciones tienen 0 observaciones, y el resultado es una matriz con ceros que la biología sabe que no son cero. El enfoque bayesiano resuelve esto sin “inventar” datos: usa explícitamente lo que ya sabes sobre el ciclo de vida (la previa) y produce estimaciones realistas con intervalos de credibilidad. Aprender a usarlo es una de las inversiones de mayor retorno en la dinámica poblacional aplicada.
Este capítulo presenta una alternativa bayesiana para estimar los parámetros de las matrices de transición y fecundidad en modelos poblacionales, especialmente útil en el estudio de especies raras. Estas especies suelen tener tamaños de muestra reducidos o transiciones vitales poco frecuentes, lo que dificulta obtener estimaciones confiables. El enfoque bayesiano propuesto permite incorporar conocimiento previo sobre el ciclo de vida de la especie para generar estimaciones más realistas.
Consideremos una especie hipotética con tres etapas en su ciclo de vida: semillas, plántulas y adultos. En este caso, sólo los adultos producen semillas. El primer paso del análisis consiste en identificar qué etapas tienen mayor influencia en la dinámica poblacional y cómo optimizar la recolección de datos en campo.
La primera figura muestra el ciclo de vida considerado. Tras definirlo, se recolectaron datos de campo y se estimaron las tasas de transición y fecundidad. Los resultados revelaron lo siguiente:
Estas interpretaciones extremas —que ningún individuo muere o que todos mueren— no son biológicamente realistas. Sin embargo, pueden surgir como resultado de un tamaño de muestra reducido. Por ejemplo, si sólo se observaron 4 plantas adultas (debido a que la especie es rara) y todas sobrevivieron, se estimaría una tasa de supervivencia del 100%, lo que implicaría que los individuos son “inmortales”. Este resultado no refleja la realidad biológica, sino una limitación del muestreo: por simple azar, no se observó ninguna muerte. Por lo tanto, es fundamental reconocer cuándo los parámetros estimados no representan adecuadamente la historia de vida típica de la especie.
raretransEl paquete raretrans fue desarrollado para abordar este tipo de situaciones. Utiliza un enfoque bayesiano que permite incorporar conocimiento previo sobre el ciclo de vida de la especie, ayudando a estimar parámetros más realistas incluso cuando los datos son escasos o incompletos.
Algunas interpretaciones biológicamente poco plausibles que pueden surgir del análisis de datos limitados incluyen:
Nota: Aunque estos resultados pueden surgir directamente de los datos recolectados, muchas veces son consecuencia de tamaños de muestra pequeños o de eventos raros que no se observaron durante el periodo de muestreo. Por ello, es esencial considerar tanto la biología de la especie como las limitaciones del muestreo al interpretar los datos.
matA <- rbind(
c(0.0, 0.0, 0.0),
c(1.0, 0.0, 0.0),
c(0.0, 0.0, 0.9)
)
stages <- c("semillas", "plantulas", "adultos")
title <- NULLEl ciclo de vida propuesto para esta especie hipotética presenta limitaciones importantes que lo hacen poco realista. En particular:
Como resultado, la tasa de crecimiento poblacional asintótica observada es \(\lambda = 0.9\) lo que implica una población en declive.
Además, la matriz poblacional estimada presenta dos propiedades estructurales problemáticas:
No es ergódica: No es posible transitar desde cualquier estado a cualquier otro estado dentro del ciclo de vida. Esto limita la conectividad entre etapas y reduce la capacidad de la matriz para representar adecuadamente la dinámica poblacional.
Es reducible: Una o más filas y columnas pueden eliminarse sin alterar las propiedades fundamentales de la matriz (como sus eigenvalores y eigenvectores). Esto indica que ciertas etapas del ciclo de vida no contribuyen significativamente a la dinámica poblacional, lo cual puede ser un artefacto del muestreo o de una estructura biológica mal representada.
plot_life_cycle(matA, stages = stages)Al estudiar especies longevas como Sequoia sempervirens, es común enfrentar interpretaciones erróneas derivadas de datos limitados o mal contextualizados. Esta especie puede alcanzar tamaños extraordinarios (diámetros a la altura del pecho, dbh, mayores a 3 metros) y tener una longevidad estimada entre 1,200 y 2,200 años.
Supongamos que se realiza un muestreo de más de 1,000 individuos grandes. Es posible que, en un intervalo de uno o incluso diez años, no se observe la muerte de ningún individuo. A partir de estos datos, se podría concluir erróneamente que la supervivencia anual es del 100%. Sin embargo, esta interpretación es biológicamente poco realista por varias razones:
Tamaño de muestra insuficiente: Si el muestreo se ampliara a 10,000 o 100,000 individuos, probablemente se detectarían muertes, aunque fueran pocas. Esto sugiere que la mortalidad es baja, pero no nula.
Relación entre tamaño y mortalidad: En plantas, la probabilidad de morir suele disminuir con el tamaño del individuo. Los árboles más grandes tienen mayor resistencia a factores ambientales y biológicos, lo que reduce su tasa de mortalidad.
Dependencia temporal de la mortalidad: La muerte de individuos puede estar influenciada por eventos bióticos (plagas, enfermedades) o abióticos (sequías, incendios, tormentas) que varían entre años. En años extremos, la mortalidad puede aumentar significativamente.
Por lo tanto, aunque los datos observados sugieren una supervivencia perfecta, la probabilidad real de mortalidad es baja pero no cero. Este tipo de sesgo en la interpretación puede ser corregido mediante enfoques estadísticos más robustos, como el uso de modelos bayesianos que incorporan incertidumbre y conocimiento previo sobre la biología de la especie.
raretransEn la siguiente sección usaremos el paquete raretrans, que facilita los estimados de las entradas en las matrices de transición y fecundidades. El manuscrito (Tremblay et al., 2021; Tyre, Tremblay & Perez, 2024) contiene los datos originales y el uso del paquete. Para tener acceso a los documentos sigan el siguiente enlace: https://doi.org/10.1016/j.ecolmodel.2021.109526.
Population projections from holey matrices: Using prior information to estimate rare transition events
Resumen
Las matrices de proyección poblacional (MPP) son una herramienta común para predecir el comportamiento de la población a corto y largo plazo. Sin embargo, cuando se utilizan para especies raras, amenazadas y en peligro de extinción los datos disponibles suelen tener tamaños de muestra reducido. En consecuencia, algunos eventos demográficos podrían perderse en la modelización dando como resultado matrices con eventos faltantes que reflejen trayectorias de ciclo de vida incompletas o biológicamente poco realistas. Estas matrices con eventos faltantes (ceros; e.g., sin observación de semillas que se transforman en plántulas) a menudo se sustituyen utilizando información de la literatura aunque información referente a otras poblaciones, otros períodos de tiempo, otras especies usualmente se omite.
Una forma de abordar este problema es a través de métodos bayesianos para el ajuste de los parámetros matriciales. Este marco estadístico nos permite usar información previa para completar los valores faltantes usando, por ejemplo, una distribución de Dirichlet para incluir información sobre las transiciones de estado y una distribución Gamma como previa para la fecundidad. Este método permite integrar información conocida (a priori) e incluir explícitamente el peso de la información previa en las distribuciones posteriores. Además, al establecer explícitamente las distribuciones previas, los resultados son reproducibles y se pueden volver a evaluar con diferentes distribuciones a priori en el futuro. A continuación mostramos cómo el uso de estos métodos en datos reales influye en la dispersión de los resultados posteriores y origina matrices irreducibles y ergódicas, permitiendo hacer inferencias biológicamente más realistas sobre las probabilidades de transición y fecundidad.
El paquete raretrans ya está disponible en CRAN, por lo que se puede instalar de la manera habitual. Solo es necesario ejecutar la instrucción una vez; en ejecuciones posteriores basta con cargar la biblioteca con library(raretrans).
# Instalar desde CRAN (ejecutar solo una vez)
install.packages("raretrans")library(raretrans)Ver el siguiente enlace para más información, la información que se presenta a continuación es una traducción y ampliación de la información del siguiente enlace. A continuación se muestra el uso del paquete raretrans para estimar los parámetros en una población y en un periodo de transición determinados.
https://atyre2.github.io/raretrans/articles/onepopperiod.html
A continuación se ejemplifica el uso del paquete raretrans para estimar los parámetros en una población y en un periodo de transición determinados.
# Mi tema de formato de las gráficas de ggplot2 personal
rlt_theme <- theme(axis.title.y = element_text(colour="grey20",size=15,face="bold"),
axis.text.x = element_text(colour="grey20",size=15, face="bold"),
axis.text.y = element_text(colour="grey20",size=15,face="bold"),
axis.title.x = element_text(colour="grey20",size=15,face="bold"))raretrans supone que la matriz de proyección es una lista con dos componentes:
Matriz de transición: Probabilidades de pasar de una etapa a otra (o permanecer en la misma).
Matriz de fecundidad: Número de nuevos individuos producidos por cada etapa reproductiva.
Este es el formato de salida de popbio::projection.matrix. Si tenemos transiciones individuales en nuestra base de datos (con clase data.frame()), podemos usar raretrans::get_state_vector para obtener el número inicial de individuos por etapa.
Podemos utilizar la función popbio::projection.matrix para obtener los datos necesarios. Hacemos una demostración con los datos de transición y fecundidad de Lepanthes eltoroensis POPNUM 250 en el periodo 5. Lepanthes eltoroensis es una orquídea epífita, endémica de Puerto Rico y con una distribución limitada a una pequeña región de la isla.
raretransdata("L_elto") # carga el conjunto de datos `L_elto` en la memoria de la computadora (los datos están incluido en el paquete `raretrans`)
head(L_elto) |> flextable()Error in `flextable()`:
! could not find function "flextable"Supongamos que tienes una base de datos llamada L_elto con las siguientes columnas:
pop: identificador de la poblaciónperiod: periodo de tiempostage: etapa actual del individuonext_stage: etapa siguientefertility: indicador de fecundidadSeleccionamos la población POPNUM 250 y el periodo 5:
El primer paso es seleccionar los datos de la población y el periodo de tiempo que queremos analizar. Posteriormente, cambiamos el nombre del estado más pequeño a “seedling”. Este paso no es estrictamente necesario pero nos ayuda a entender más fácilmente la nomenclatura.
La base de datos L_elto está compuesta por datos individuales del estado actual (stage, periodo t), del estado siguiente (next_stage, periodo t+1) y de la fertilidad. Tenga en cuenta que aquí “p” significa “plántula” por la traducción al español. En la siguiente caja de código las primeras líneas seleccionan la población de interés y el periodo de tiempo sobre el que queremos trabajar y cambian el nombre de la etapa del ciclo vital de “p” a “s”.
onepop <- L_elto %>%
filter(POPNUM == 250, year == 5) %>% # Filtrar la población # 250, el periodo (año=year) 5
mutate(
stage = case_when(stage == "p" ~ "s", TRUE ~ stage),
next_stage = case_when(next_stage == "p" ~ "s", TRUE ~ next_stage)
) # redefine "p" por plantula a "s" para seedlingEn el siguiente código usamos la función popbio::projection.matrix() para crear la matriz de proyección poblacional. Tenemos ahora nuestra matriz de transición y fecundidad basadas en los datos de la población #250 del periodo 5. La función projection.matrix() requiere que se especifiquen los nombres de las etapas y los estados siguientes, así como el nombre de la columna de fecundidad. En este caso, las etapas son “s” (plántula), “j” (juvenil) y “a” (adulto), y el estado siguiente es next_stage. La columna de fecundidad se llama fertility. El argumento sort especifica el orden en que queremos que aparezcan las etapas en la matriz.
Solamente los adultos pueden producir plántulas, por lo que la matriz de fecundidad tiene un valor solamente en la posición (1,3), (plántula, adulto). La función projection.matrix() devuelve una lista con dos matrices: una de transición y otra de fecundidad. La matriz de transición es una matriz cuadrada que muestra las probabilidades de transición entre las etapas, mientras que la matriz de fecundidad es una matriz que muestra la fecundidad por etapa.
A partir de este punto los análisis se llevan a cabo considerando s (plántula), j (juvenil) y a (adulto), y las dos matrices: T (transición de estadios) y F (fecundidad). La tasa de crecimiento asintótica de la población observada es \(\lambda\). Las transiciones raras que faltan en nuestra primera matriz de transición TF$T, pero que sabemos que ocurren, son la transición de plántula (s) a adulto (a) y la transición de j a s.
# popbio::projection.matrix no funciona con el formato *tibbles*, por consecuencia se convierte en data.frame
head(onepop) |> flextable() # Ahora tenemos solamente datos de la población #250 del periodo 5Error in `flextable()`:
! could not find function "flextable"# Crear TF = TRUE, añadir para formatear correctamente.
TF <- popbio::projection.matrix(
as.data.frame(onepop),
stage = stage,
fate = next_stage,
fertility = "fertility",
sort = c("s", "j", "a"),
TF = TRUE
)
TF # Esta es la estructura de etapas de vida para esa población. Nota que tenemos dos matrices, una de transiciones **T** y otra de fecundidad **F**.$T
s j a
s 0.09090909 0.00000000 0.00000000
j 0.63636364 0.57446809 0.00000000
a 0.00000000 0.29787234 0.85294118
$F
s j a
s 0.0000000 0.0000000 0.1176471
j 0.0000000 0.0000000 0.0000000
a 0.0000000 0.0000000 0.0000000Dado que nuestros conteos (número de individuos), \(N\) y el tamaño de muestra equivalente a priori se expresan como múltiplos del número de individuos observados, necesitamos obtener el número de individuos en cada etapa (\(N\)) en el primer periodo de tiempo. Para esto utilizamos la función raretrans::get_state_vector(). Se observa que la cantidad de individuos iniciales por etapa es el número de individuos en el primer muestreo es muy limitada, con 11 plántulas, 47 juveniles y 34 adultos, para esta población en este periodo de muestreo.
N <- get_state_vector(onepop, stage = stage, sort = c("s", "j", "a"))
N # Un vector # de individuos iniciales para cada etapa, nota que la cantidad por etapa, "stage", son los individuos en el primer muestreo[1] 11 47 34La lista de matrices y el vector de conteos de individuos no tienen por qué proceder de un data frame como hemos hecho anteriormente. Mientras tengan el formato esperado, pueden crearse a mano. Usemos la población 231 en el periodo 2 como ejemplo para dividir la matriz poblacional en las submatrices de transición T y de fecundidad F. En el siguiente código, m denota mortalidad, es decir, plantas que están muertas. Note que aquí las transiciones están estimadas con un tamaño de muestra todavía más pequeño: solamente 2 plántulas, 6 juveniles y 16 adultos. Ninguna de las 2 plántulas pasó a la etapa de juveniles; en consecuencia, esta transición es cero. Además, una de estas plántulas murió, por lo que la supervivencia es del 50%. Sin embargo, si hubieran muerto las dos plántulas la mortalidad sería del 100%, y si hubieran sobrevivido las dos plántulas la supervivencia sería del 100%: ¡no hay valores intermedios!
La matriz de transición, la matriz de fecundidad y el vector de tamaño poblacional para la población 231, periodo 2, son:
TF2$Tmat
stage
fate p j a
p 0.5000000 0.0000000 0.0000000
j 0.0000000 0.8333333 0.0000000
a 0.0000000 0.1666667 0.8750000
$Fmat
[,1] [,2] [,3]
[1,] 0 0 0.125
[2,] 0 0 0.000
[3,] 0 0 0.000N2 p j a
2 6 16 A la matriz mostrada le falta la transición de plántula a juvenil. Además, dado que ninguno de los 6 juveniles murió, hay una sobreestimación de su supervivencia. La tasa de crecimiento asintótica de la población observada es \(\lambda\). La matriz no es ergódica (no se puede llegar a cualquier otro estado desde uno o más estados), y es reducible, lo que significa que una o más columnas y filas se pueden descartar y tienen las mismas propiedades en términos de sus eigenvalores y eigenvectores.
Lo que vimos anteriormente fueron algunos de los problemas reales que encontramos cuando usamos datos de campo y hay tamaños de muestra pequeños o transiciones raras. A continuación veremos cómo resolverlos.
stages2 <- c("plantulas", "juveniles", "adultos")
title <- NULL
plot_life_cycle(Tmat, stages = stages2)Para evitar que el modelo genere transiciones imposibles entre etapas del ciclo de vida, se pueden incorporar previas informativas. Estas previas se basan en el conocimiento experto sobre la especie, en este caso proporcionado por un especialista en orquídeas epífitas del género Lepanthes (RLT). La información debe tener el mismo formato que la matriz de transición, es decir:
Es importante que estos valores se definan antes de recolectar los datos de campo, para evitar que el análisis se vea influenciado o sesgado por los datos observados.
RLT_Tprior <- matrix(c(0.25, 0.025, 0.0,
0.05, 0.9, 0.025,
0.01, 0.025, 0.95,
0.69, 0.05, 0.025),
byrow = TRUE, nrow = 4, ncol = 3)Note que en la matriz, la 1ª fila, 3ª columna es 0.0, porque esta transición es imposible. Las previas se construyen de manera que las columnas sumen 1, lo que crea mayor flexibilidad para la ponderación de los valores a priori. Por defecto, la suma es 1, interpretado como un tamaño de muestra a priori de 1.
Ahora, usando la función fill_transitions se agregan los valores del campo TF, el tamaño de muestra en el tiempo inicial N, la matriz de valores con priors informativos, y se calcula la matriz de transición a posteriori. Compare TF con esta matriz para ver los cambios.
fill_transitions(TF, N, P = RLT_Tprior) [,1] [,2] [,3]
[1,] 0.1041666667 0.0005208333 0.0000000000
[2,] 0.5875000000 0.5812500000 0.0007142857
[3,] 0.0008333333 0.2921875000 0.8557142857Cuando usamos previas informativas en el modelo, podemos ajustar su peso para reflejar cuánto creemos en la previa frente a los datos. Esto se hace con el argumento priorweight de fill_transitions().
¿Cómo funciona priorweight exactamente? El argumento es un multiplicador sobre el tamaño de muestra observado \(N\):
priorweight × \(N\) (en unidades de “observaciones equivalentes”).priorweight × \(N\) = \(N \times (1 + \text{priorweight})\).priorweight = -1, lo que hace que la previa contribuya solamente su suma original (típicamente 1, ya que cada columna de la matriz de previas se construye para sumar 1).fill_transitions(TF, N, P = RLT_Tprior, priorweight = 0.5) [,1] [,2] [,3]
[1,] 0.143939394 0.008333333 0.000000000
[2,] 0.440909091 0.682978723 0.008333333
[3,] 0.003333333 0.206914894 0.885294118El “peso de la previa” es un tamaño de muestra equivalente. Asignar priorweight = 0.5 con \(N = 2\) observaciones es como decir: “antes de salir al campo, mi creencia tiene la misma fuerza que \(0.5 \times 2 = 1\) observación adicional”. Si tu previa tiene peso efectivo 5 y observas 50 transiciones, los datos pesarán 10 veces más que la previa. La regla práctica es dejar que los datos dominen: el peso de la previa debería ser mucho menor que el tamaño de muestra observado, salvo que tengas razones biológicas muy fuertes para confiar en la previa.
Algunos valores típicos de priorweight:
priorweight = -1 (defecto): la previa contribuye su suma original, típicamente 1 observación equivalente. Útil como previa débil.priorweight = 0.1 a 0.5: la previa pesa entre 10% y 50% del tamaño de muestra observado. Apropiado cuando hay conocimiento experto pero los datos deben dominar.priorweight = 1: la previa pesa lo mismo que los datos. Sólo justificado si los datos son muy escasos y la previa es muy confiable.priorweight > 1: la previa pesa más que los datos. Rara vez justificado en estudios poblacionales.Cuándo no usar una previa informativa fuerte. Las previas informativas son apropiadas cuando hay conocimiento biológico sólido que no aparece en los datos (e.g., experiencia directa con la especie, estudios en poblaciones cercanas). No son apropiadas cuando:
En estos casos, conviene usar previas débiles (poco peso) o reportar resultados con y sin previa para mostrar la sensibilidad a la elección.
Recomendación general: asigna a las previas un peso menor (preferiblemente mucho menor) que el tamaño de muestra observado. Esto garantiza que los datos reales tengan mayor peso en el análisis y evita que el conocimiento previo domine injustamente los resultados.
Resultado
La nueva matriz de transición:
Una vez ajustada la matriz de transición con los datos observados y las previas informativas, podemos calcular los intervalos de credibilidad (ICr) para cada entrada de la matriz.
Por qué Dirichlet (transiciones) y Gamma (fecundidad). Cada columna de una matriz de transición es un vector de probabilidades que suman 1. La distribución natural sobre tales vectores es la distribución Dirichlet, generalización multidimensional de la Beta. Al observar transiciones (datos multinomiales) y combinar con una previa Dirichlet, la posterior también es Dirichlet (es una previa conjugada) — esto hace los cálculos exactos y rápidos. Para cada entrada individual de la matriz, la distribución marginal es Beta, lo cual permite calcular intervalos de credibilidad directamente con qbeta(). La fecundidad, por ser un conteo positivo no acotado, se modela con una previa Gamma sobre la tasa de Poisson.
Intervalo de confianza (IC) versus intervalo de credibilidad (ICr). Aunque suenan parecidos, significan cosas distintas:
Lo que la mayoría de las personas intuitivamente quieren decir al leer un IC frecuentista es la interpretación bayesiana — y los intervalos de credibilidad son los que efectivamente se la dan.
¿Cómo se calculan? Cada elemento de la matriz de transición ajustada sigue una distribución beta, ya que proviene de una distribución posterior marginal de una Dirichlet (multinomial-Dirichlet). Esto permite calcular intervalos de credibilidad para las tasas de transición.
Para obtener estos valores, se puede usar el tipo de retorno TN, que proporciona:
Supongamos que estamos analizando las transiciones desde la etapa de plántula (s). Si el número de transiciones observadas es 2, y el peso de las previas es 1, el tamaño de muestra efectivo es 3. Esto genera un intervalo de credibilidad más amplio, reflejando la alta incertidumbre en la estimación. Por ejemplo:
Este rango amplio indica que el mejor estimado es muy incierto, lo cual es típico cuando se trabaja con tamaños de muestra pequeños.
Un intervalo de credibilidad amplio significa que hay mucha incertidumbre en la estimación del parámetro. Esto resalta la importancia de contar con:
TN <- fill_transitions(
TF,
N,
P = RLT_Tprior,
priorweight = 0.1,
returnType = "TN"
)
a <- TN[, 1] # cambie 1 a 2, 3 etc para obtener la distribución beta marginal de cada columna (etapa de vida).
b <- sum(TN[, 1]) - TN[, 1] # cambie 1 a 2, 3 etc para obtener la distribución beta marginal de cada columna.
p <- a / (a + b)
lcl <- qbeta(0.025, a, b) # el valor 0.025 de la distribución beta
ucl <- qbeta(0.975, a, b) # el valor 0.975 de la distribución beta
knitr::kable(sprintf("%.4f (%.4f, %.4f)", p, lcl, ucl)) # Crear una tabla con los resultados con 4 cifras decimales (note la ultima fila son los fallecidos)| x |
|---|
| 0.1054 (0.0058, 0.3212) |
| 0.5831 (0.3087, 0.8315) |
| 0.0009 (0.0000, 0.0052) |
| 0.3107 (0.0956, 0.5846) |
Una de las grandes ventajas de usar el paquete raretrans es que podemos calcular los intervalos de credibilidad de las tasas de transición y fecundidad de manera rápida y sencilla. Además, los intervalos de credibilidad están limitados a los valores entre 0 y 1, lo que es importante para las tasas de transición y fecundidad. Esto es especialmente útil cuando se trabaja con datos de campo donde el tamaño de muestra es pequeño y las transiciones raras son comunes. Muchos métodos estadísticos no pueden calcular intervalos de confianza para tasas de transición y fecundidad, y cuando lo hacen tienen valores ilógicos, como proporciones negativas o mayores a 1.
Cuando se aumenta el tamaño de muestra efectivo, los intervalos de credibilidad se vuelven más estrechos, lo que indica mayor confianza en las estimaciones puntuales. Aumentar el tamaño de muestra efectivo a \(20\) especificando: priorweight\(= 9 (9*2 = 18 + 2 = 20)\) hace que los intervalos de credibilidad se reduzcan bastante. Entonces, es importante enfatizar en que el tamaño de muestra tiene gran impacto sobre la confianza que se tiene sobre el estimador puntual (el promedio) de las transiciones, permanencias y mortalidad. Note en la figura como el tamaño del intervalo de credibilidad cambia drásticamente cuando se aumenta el número de datos.
El intervalo de credibilidad se estrecha cuando el tamaño de muestra efectivo aumenta, reflejando mayor confianza en el estimador puntual.
IMPORTANTE: Aunque este ejemplo muestra cómo el tamaño de muestra afecta la credibilidad, en la práctica:
TN <- fill_transitions(
TF,
N,
P = RLT_Tprior,
priorweight = 9,
returnType = "TN"
)
a <- TN[, 1]
b <- sum(TN[, 1]) - TN[, 1] # El parámetro beta es la suma de todos los demás parámetros de Dirichlet.
p2 <- a / (a + b)
lcl2 <- qbeta(0.025, a, b)
ucl2 <- qbeta(0.975, a, b)
alltogethernow <- tibble(
priorweight = factor(rep(c(1, 9), each = 4)),
fate = factor(
rep(c("s", "j", "a", "dead"), times = 2),
levels = c("s", "j", "a", "dead")
),
p = c(p, p2),
lcl = c(lcl, lcl2),
ucl = c(ucl, ucl2)
)
ggplot(
data = alltogethernow,
mapping = aes(x = fate, y = p, group = priorweight, colour = priorweight)
) +
geom_point(
mapping = aes(shape = priorweight),
position = position_dodge(width = 0.5)
) +
geom_errorbar(
mapping = aes(ymin = lcl, ymax = ucl),
position = position_dodge(width = 0.5)
) +
theme_bw() +
rlt_style_colour()
Obtener intervalos de credibilidad para la tasa de crecimiento asintótica, \(\lambda\), requiere simular matrices a partir de las distribuciones posteriores. Estas simulaciones usualmente son complicadas de hacer, por lo que hemos incorporado en el paquete la función raretrans::sim_transitions() para generar una lista de matrices simuladas dada la matriz observada y las previas. En otras palabras, esta función calcula múltiples veces las matrices de transición y fecundidad basándose en los datos observados y en los valores a priori especificados, pero variando los valores de los parámetros de las distribuciones beta y gamma.
Siguiendo el marco bayesiano que hemos venido trabajando en este capítulo, es necesario especificar las previas de la fecundidad como una matriz. Las celdas de la matriz de los estados en los que no hay reproducción deben ser NA, usando NA_real_, el cual es un valor que no está presente pero podría aceptar valores con puntos decimales. Nota que el valor a priori de la fecundidad es 0.0001.
alpha <- matrix(c(NA_real_, NA_real_, 1e-5,
NA_real_, NA_real_, NA_real_,
NA_real_, NA_real_, NA_real_), nrow=3, ncol = 3, byrow = TRUE)
beta <- matrix(c(NA_real_, NA_real_, 1e-5,
NA_real_, NA_real_, NA_real_,
NA_real_, NA_real_, NA_real_), nrow=3, ncol = 3, byrow = TRUE)
fill_fertility(TF, N, alpha = alpha, beta = beta)
s j a
s 0.0000000 0.0000000 0.1176473
j 0.0000000 0.0000000 0.0000000
a 0.0000000 0.0000000 0.0000000El cambio en la fecundidad es < 0.0001 en comparación con el valor observado.
Corra el script múltiples veces, y verá que los valores cambian en cada simulación. En este script la simulación es solamente una vez (samples = 1). Este script es para que visualicen lo que hace la función sim_transitions().
sim_transitions(
TF,
N,
P = RLT_Tprior,
alpha = alpha,
beta = beta,
priorweight = 0.5,
samples = 1
)[[1]]
[,1] [,2] [,3]
[1,] 0.05605405 0.005638743 0.218884826
[2,] 0.31517061 0.616246548 0.006185498
[3,] 0.02258141 0.221289566 0.918125802Ahora simulamos 5000 veces, calculamos el valor \(\lambda\) de cada matriz y creamos un histograma de su distribución. Se añaden los valores de lambda en un objeto llamado RLT_0.5. Se visualiza la distribución de \(\lambda\) usando un histograma. Note que el valor de \(\lambda = 1\) queda dentro de la distribución; en otras palabras, el histograma muestra que la mayoría de los valores de \(\lambda\) están cerca de 1.
#set.seed(8390278) # make this part reproducible
alpha2 <- matrix(c(NA_real_, NA_real_, 0.025,
NA_real_, NA_real_, NA_real_,
NA_real_, NA_real_, NA_real_), nrow=3, ncol = 3, byrow = TRUE)
beta2 <- matrix(c(NA_real_, NA_real_, 1,
NA_real_, NA_real_, NA_real_,
NA_real_, NA_real_, NA_real_), nrow=3, ncol = 3, byrow = TRUE)
# generar 5000 matrices basadas en las previas de transiciones y de fecundidades, el tamaño de muestra, además de los datos
RLT_0.5s <- sim_transitions(TF, N, P = RLT_Tprior, alpha = alpha2, beta = beta2,
priorweight = 0.5, samples = 5000)
# extraer los valores de lambda de cada simulación en un histograma
RLT_0.5 <- tibble(lposterior = map_dbl(RLT_0.5s, lambda)) # convertir la lista en un tibble
ggplot(data = RLT_0.5,
mapping = aes(x = lposterior)) +
geom_histogram(binwidth = 0.01, colour="white") +
theme_bw() +
rlt_style_fill() +
geom_vline(xintercept = 1,
color = "#009392", linewidth=1.5)
A diferencia de un contraste de hipótesis frecuentista (donde se obtiene un valor-p), el marco bayesiano nos permite responder directamente la pregunta de interés: ¿con qué probabilidad la población está creciendo? Lo hacemos contando qué proporción de las simulaciones de la distribución posterior produjo \(\lambda > 1\). En este caso, \(P(\lambda > 1 \mid \text{datos}, \text{previas}) = 0.14\), lo que indica que sólo en el 14% de los escenarios consistentes con los datos y la previa la población crece. La evidencia para crecimiento poblacional es débil; lo más probable es que la población esté en declive. La estadística p_increase reporta exactamente este valor.
RLT_0.5_summary <- summarize(
RLT_0.5,
medianaL = median(lposterior),
promedioL = mean(lposterior),
skewnessL=e1071::skewness(lposterior, type=2),
ICrBajo = quantile(lposterior, probs = 0.025),
ICrAlto = quantile(lposterior, probs = 0.975),
p_increase = sum(lposterior > 1.) / n()
)flextable(signif(RLT_0.5_summary, digits = 2))Error in `flextable()`:
! could not find function "flextable"Esta sección es para demostrar por qué no es recomendable el uso de distribuciones previas uniformes. Si se entiende el concepto de valores a priori y a posteriori, no es necesario leer esta sección.
A continuación vamos a añadir una distribución Dirichlet uniforme como previa con un peso = \(1\) a la matriz de transición, \(T\). En este caso, tenemos 4 destinos (3 + muerte), por lo que cada destino añade 0.25 a la matriz de destinos observados (¡no a la matriz de transición!). Cuando especificamos una matriz con valores a priori para las transiciones, hay una fila más que columnas. Esta fila extra representa la mortalidad por etapa.
Note que no se recomienda el uso de una previa uniforme; aquí se utiliza sólo para mostrar el concepto. Una distribución previa uniforme no refleja un ciclo de vida razonable para especies biológicas. Por ejemplo, una previa uniforme podría hacer que las plántulas pasen a adulto de un tiempo a otro o que los adultos regresen a ser plántulas, lo cual no tiene mucho sentido biológico en muchas especies. Considera lo siguiente: ¿es posible que una ballena bebé pueda pasar a ser madre al año siguiente? Añadir una transición de 0.25 de ballena bebé a madre no tiene sentido biológico.
Note que ahora la matriz de transición es ergódica pero no tiene sentido biológico. Compare la matriz con las previas uniformes versus la matriz original.
Tprior <- matrix(0.25, byrow = TRUE, ncol = 3, nrow = 4) # crear la matriz de previas uniformes
Tprior [,1] [,2] [,3]
[1,] 0.25 0.25 0.25
[2,] 0.25 0.25 0.25
[3,] 0.25 0.25 0.25
[4,] 0.25 0.25 0.25Posterior_con_Unif = fill_transitions(TF, N, P = Tprior) # resultado de la matriz de transición básica con las previas
Posterior_con_Unif [,1] [,2] [,3]
[1,] 0.10416667 0.005208333 0.007142857
[2,] 0.60416667 0.567708333 0.007142857
[3,] 0.02083333 0.296875000 0.835714286# Para entender las diferencias del impacto de las matriz compara los resultados con *$T* del objeto *TF*
TF$T
s j a
s 0.09090909 0.00000000 0.00000000
j 0.63636364 0.57446809 0.00000000
a 0.00000000 0.29787234 0.85294118
$F
s j a
s 0.0000000 0.0000000 0.1176471
j 0.0000000 0.0000000 0.0000000
a 0.0000000 0.0000000 0.0000000Podemos obtener el resultado haciendo los cálculos directamente. Lo primero que necesitamos es el vector de observaciones porque la distribución posterior se calcula a partir de las observaciones de transiciones, no de la matriz de transiciones. Si no tiene experiencia con análisis Bayesiano es recomendable seguir los siguientes scripts
Tobs <- sweep(TF$T, 2, N, "*") # obtener las observaciones de transiciones
Tobs <- rbind(Tobs, N - colSums(Tobs)) # añadir la fila de muerte
Tobs <- Tobs + 0.25 # añadir los valores de las previas
sweep(Tobs, 2, colSums(Tobs), "/")[-4, ] # dividir por la suma de la columna y descarta la fila de muerte s j a
s 0.10416667 0.005208333 0.007142857
j 0.60416667 0.567708333 0.007142857
a 0.02083333 0.296875000 0.835714286La previa uniforme “rellena” las transiciones faltantes, pero también crea problemas porque algunos de los valores de transición que se asignan son biológicamente imposibles en muchas especies. Por ejemplo, para la transición de adulto a plántula se asigna un valor; sin embargo, este valor sólo es posible en la matriz de fecundidad \(F\) en las especies de Lepanthes. Debido a esto, no es recomendable el uso de previas uniformes, ya que no toman en cuenta el ciclo de vida de nuestra especie. Será necesario crear una matriz de priors que refleje el ciclo de vida de la especie de interés.