Clicky

logo

Introducción

La primera parte de este tutorial tratará sobre cómo transformar un modelo frecuentista básico en un modelo bayesiano e introducirá algunos conceptos como el grafo acíclico dirigido (DAG, por sus siglas en inglés) y la distribución a priori.

La segunda parte se dedicará a introducir la modelización de poblaciones a partir de una muestra de sitios de estudio escasa pero representativa. Presentaremos el conjunto de datos con el que trabajaremos durante los tres próximos tutoriales y crearemos nuestros dos primeros modelos de población:

  • un modelo basado en la distribución normal.

  • un modelo multinivel basado en un compuesto Poisson lognormal.

Será la ocasión de experimentar con el software Stan y su interfaz R rstan.

Objetivos

  1. Escribir una regresión lineal simple en un marco bayesiano

  2. Adaptar las herramientas del estadístico a un ejemplo real

  3. Ajustar un modelo normal en stan con datos simulados:

    1. Formato de datos parastan

    2. Especificar un modelo en el lenguajestan

    3. Configurar una muestra MCMC para ajustar el modelo

  4. Ajustar un modelo de Poisson para el recuento de población

    1. Evaluar los resultados y las limitaciones

  5. Ajustar un modelo Poisson lognormal para modelizar la población

Lecturas complementarias

Estos tutoriales no son una introducción a la estadística. Para esta lección en concreto, sería bueno familiarizarse con algunos conceptos estadísticos y, para ello, indicamos recursos útiles:

Y añadimos a esa lista la documentación del software utilizado:

  • Software: Stan

    Stan es una biblioteca C++ para modelización e inferencia bayesiana que utiliza principalmente el algoritmo No-U-Turn sampler (NUTS), introducido por Hoffman y Gelman (s. f.) para obtener simulaciones posteriores con un modelo y datos especificados por el usuario

  • Interfaz: rstan

    El paquete rstan permite ajustar cómodamente modelos de Stan desde R (R Core Team 2014) y acceder a los resultados, incluidas inferencias posteriores y cantidades intermedias como evaluaciones de la densidad posterior logarítmica y sus gradientes.

De la mentalidad frecuentista a la bayesiana

En un enfoque frecuentista estándar, una regresión lineal entre \(Y\), la variable de respuesta, y \(X\), los indicadores, puede formularse como:

\[\begin{equation} Y = \alpha + \beta X + \epsilon \\ \epsilon \sim Normal(mean=0,sd=\sigma)\tag{1} \end{equation}\]

La ecuación (1) puede reescribirse como: \[Y \sim Normal(mean=\mu, sd=\sigma)\] \[ \mu=\alpha + \beta X\]

Este formato es más flexible cuando empezamos a trabajar con estructuras de error no normales y componentes de modelización personalizados.

Esta regresión lineal puede representarse mediante un grafo acíclico dirigido (DAG, por sus siglas en inglés) que ayuda a representar las relaciones entre los parámetros del modelo y los datos de entrada:

En un DAG, los cuadrados representan datos y los círculos, parámetros. Las direcciones de las flechas indican dependencia.

En un modelo bayesiano, todos los parámetros del nodo raíz (los que no tienen flechas que apunten hacia ellos) necesitan que se especifique la distribución a priori:

\[ \alpha \sim Normal(mean=0, sd=1000)\] \[ \beta \sim Normal(mean=0, sd=1000)\] \[ \sigma \sim Uniform(min=0, max=1000)\]

Estos son ejemplos de distribución a priori débilmente informativa porque las varianzas son grandes en relación con los datos. La distribución a priori poco informativa no debería tener ninguna influencia notable en las estimaciones finales de los parámetros.

Cómo elegir la distribución a priori

Para identificar una distribución a utilizar para las variables a priori, pregúntate:

  1. "¿Qué valores son posibles para este parámetro?" y

  2. "¿Qué tan seguro estoy de ello?"

En la mayoría de los análisis bayesianos, generalmente sabemos cómo son y cómo no son nuestros datos, de modo que podemos construir distribuciones a priori que reflejen las propiedades de los conjuntos de datos potenciales. La definición de tales distribuciones a priori consiste en indicar cuál es el rango probable de valores para nuestros parámetros y reducir la ponderación de los valores fuera de rango. Además, si no estamos seguros del rango, podemos aumentar la variabilidad de la distribución a priori (a través de su parámetro de escala/varianza) porque no queremos que las variables a priori influyan demasiado en las variables a posteriori.

Los coeficientes de regresión suelen ser números continuos que pueden tomar valores de -∞ a +∞. La distribución normal es una buena elección a priori para estos parámetros porque tiene las mismas características.

Las desviaciones típicas son números continuos que deben ser positivos. Una distribución normal no es una buena elección a priori en este caso porque incluye números negativos. La más sencilla es una distribución uniforme (0, a) con a definiendo el límite superior de los valores probables. En realidad, esto no es perfecto porque impone una restricción dura. Otras distribuciones pueden considerarse como una distribución media normal o media de Cauchy.

El libro Bayesrules! tiene un muy buen capítulo sobre la interacción entre la distribución a priori y los datos.

El equipo de stan ha elaborado interesantes directrices para el ajuste previo.

Probemos con datos simulados

Configuración

En primer lugar, descarga las versiones más recientes de los siguientes programas:

A continuación, instala y configura rstan siguiendo atentamente las instrucciones.

# stan setup
options(mc.cores = parallel::detectCores()) #set up the maximum number of cores used by stan
rstan::rstan_options(auto_write = TRUE) # speed up running time
  • Instala un conjunto de paquetes de gestión de datos:
install.packages(c( "tidyverse", # for data manipulation
                    'kableExtra', # for good visualisation of tables
                    'here' # for handling relative path to external assets (pic, data)
                    ), 
                 dependencies = TRUE)

Para más información, consulta sus viñetas: tidyverse, kableExtra y aquí.

Simulación de datos

Simularemos observaciones falsas para introducir los conceptos básicos de la modelización bayesiana y su aplicación en stan.

Producimos nuestros datos falsos como 1000 extracciones de una distribución normal con media 5 y desviación típica 50:

# 2 Simulated data ----
# Simulate data
seed <- 2004
set.seed(seed)
data <- tibble(y=rnorm(1e3, mean= 5, sd=50))

Obsérvese que definimos seed un patrón para que los resultados se repliquen exactamente.

La distribución observada de nuestros datos simulados es la siguiente:

# plot simulated data
ggplot(data, aes(x=y))+
  geom_histogram(bins = 50)+
  theme_minimal()+
  geom_vline(xintercept = mean(data$y), color='orange', size=1)+
  annotate('text', x=15, y=30, label=paste0('Observed mean: ', round(mean(data$y),2)), colour='orange', angle=90)+
  geom_segment(x=0,y=1, xend=sd(data$y), yend=1, colour='orange', size=1)+
  annotate('text', x=50, y=4, label=paste0('Observed sd: ', round(sd(data$y),2)), colour='orange')+
  labs(y='', x='observed y')+
  theme(axis.text.y = element_blank())
Simulated observations distribution

Imagen 1: Distribución de las observaciones simuladas

Modelización de los datos

Imaginemos que no sabemos cómo se han creado las observaciones \(y\). Queremos crear un modelo de \(y\), y vemos que tiene una clara forma de campana. Por lo tanto, se puede calcular con una distribución normal:\[Y \sim Normal(mean=\mu, sd=\sigma)\]

El DAG correspondiente es:

# Normal model DAG
dagify(
  Y ~ mu,
  Y ~ sigma) %>%
  tidy_dagitty(seed=41) %>% 
  mutate(color=c('parameter', 'parameter','data')) %>% 
  ggplot(aes(x = x, y = y, xend = xend, yend = yend, color=color,shape=color)) +
  geom_dag_point() +
  geom_dag_edges() +
  geom_dag_text(col = "grey20",size=6,  parse=T) +
  scale_shape_manual(values=c(15,19))+
  theme_dag()+ labs(title = 'Model with simulated data: Normal distribution of Y', color='', shape='')

Necesitamos dar alguna prioridad a los dos parámetros: \(\mu\) la media y span class="math inline">\(\sigma\) la desviación estándar. Imaginemos que no tenemos acceso a más información sobre la variable \(Y\) que su manifestación en nuestro conjunto de datos.

Podemos suponer que \(\mu\) ronda el cero, pero sin gran confianza (véase la imagen 1). En otras palabras, podemos suponer que \(\mu\) puede extraerse de una normal (0, 100).

# define mu prior
data$prior_mu <- rnorm(1e3, mean= 0, sd=100)

ggplot(data)+
  geom_histogram(aes(x=y, after_stat(density)), bins=50, fill='grey30')+
  geom_density(aes(x=prior_mu), colour='#00BFC4', size=1)+
  theme_minimal()+
  annotate('text', y=0.0020, x=250, label=paste0('Normal(0,100)'), colour='#00BFC4', size=5)+
    theme(axis.text.y = element_blank())+
  labs(y='', x='observed y')
Mu prior distribution

Imagen 2: Distribución a priori Mu

La imagen 2 muestra lo que significa nuestra distribución a priori para \(\mu\): pensamos que es probable que \(\mu\) ronde el cero, pero podría llegar a 200. Dado que nuestra \(y\) observada, \(\hat{y}\), no tiene ocurrencias por encima de 200, es improbable que la media de \(y\) sea 200, pero no excluimos la posibilidad. Sin embargo, consideramos que 0 es más probable que 200.

Para \(\sigma\) tenemos expectativas más estrictas. De hecho, las desviaciones típicas son positivas. Así pues, elegimos como variable a priori para \(\sigma\) una distribución uniforme (0, 200).

# define sigma prior
data$prior_sigma <- runif(1e3,min = 0, max = 200)

ggplot(data)+
  geom_histogram(aes(x=y, after_stat(density)), bins=50, fill='grey30')+
  geom_density(aes(x=prior_sigma), colour='#00BFC4', size=1, trim=T)+
  theme_minimal()+
  annotate('text', y=0.0060, x=150, label=paste0('Uniform(0,200)'), colour='#00BFC4', size=5)+
      theme(axis.text.y = element_blank())+
  labs(y='', x='observed y')
Sigma prior distribution

Imagen 3: Distribución a priori Sigma

La imagen 3 muestra que una distribución a priori uniforme implica una probabilidad igual de que \(\sigma\) sea 0 o 200.

El modelo final es:

\[\begin{equation} Y \sim Normal( \mu, \sigma ) \\[1pt] \\ \mu \sim Normal( 0, 100 ) \\ \sigma \sim Uniform( 0, 200 )\tag{2} \end{equation}\]

Aplicación del modelo en stan

Para estimar \(\mu\) y \(\sigma\), escribiremos nuestro primer modelo en stan.

// Model for simulated data: y as normal distribution
data {
  int<lower=0> n; // number of observations
  real y[n]; // observations
}

// The parameters accepted by the model. Our model
// accepts two parameters 'mu' and 'sigma'.
parameters {
  real mu;
  real<lower=0> sigma;
}

// The model to be estimated. We model the output
// 'y' to be normally distributed with mean 'mu'
// and standard deviation 'sigma'.
model {
  y ~ normal(mu, sigma);
  
  mu ~ normal(0,100);
  sigma ~ uniform(0,200);
}

Un modelo stan se compone de bloques de código. Los fundamentales son:

  • el bloque de datos que describe los datos de entrada al modelo

  • el bloque de parámetros que describe los parámetros que deben estimarse

  • el bloque del modelo que describe los elementos estocásticos: (1) la interacción entre los parámetros y los datos, (2) la distribución a priori

El software stan nos obliga a declarar todas las variables, tanto parámetros como datos, con su tipo (int, real) y tamaño (indicado con [n])1. Es posible incorporar restricciones en el soporte de la variable, por ejemplo, no es posible tener una \(\sigma\) negativa (real<lower=0> sigma).

Nota: stan requiere dejar una línea en blanco al final del script.

Almacenaremos el modelo en un archivo stan llamado tutorial1_model.stan en la carpeta tutorial1.

Preparación de los datos para stan

El software Stan toma como entrada una lista de los datos observados que define las variables indicadas en el bloque de datos.

# prepare data for stan
stan_data <- list(
  y = data$y,
  n = nrow(data))

Ejecución del modelo

Configuramos los parámetros del algoritmo de la cadena de Markov.

# mcmc settings
chains <- 4
warmup <- 250
iter <- 500

El argumento de las chains especifica el número de cadenas de Markov que deben ejecutarse simultáneamente. Queremos que las cadenas de Markov reproduzcan un proceso totalmente aleatorio. Sin embargo, el diseño del algoritmo en cadena hace que cada muestra dependa de la anterior. Para recrear un escenario aleatorio, ejecutamos independientemente varias cadenas para explorar el espacio de parámetros y que, con suerte, converjan a la misma solución consistente.

El parámetro de warmup es el número de muestras al principio del proceso de estimación que descartamos de los resultados. Esto es similar a cocinar panqueques, en el sentido de que se necesita que el algoritmo caliente antes de acercarse a valores razonables.

El parámetro iter especifica el número de iteraciones, es decir, la longitud de la cadena de Markov. Cuanto más larga sea la cadena, más probable es que se estabilice en torno a la estimación correcta.

A continuación, definimos los parámetros que queremos controlar, que se almacenan durante el proceso de estimación:

# parameters to monitor
pars <- c('mu','sigma')

¡Y ya estamos listos para ejecutar el modelo!

# mcmc
fit <- rstan::stan(file = file.path('tutorial1_model.stan'), 
                   data = stan_data,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)

Comprobación de las simulaciones MCMC

Comprobamos las cadenas de Markov para ver si convergen en una solución única. Se puede visualizar con un trazado para los dos parámetros, \(\mu\) y \(\sigma\). Un trazado describe la evolución de la estimación de los parámetros a lo largo de las iteraciones de la cadena de Markov. Un buen trazado muestra la mezcla de las distintas cadenas, prueba de la convergencia hacia una única estimación.

# plot trace
stan_trace(fit, inc_warmup = T)
Model Traceplot

Imagen 4: Modelo de trazado

La imagen 4 muestra tanto el período de calentamiento (hasta 250) como las iteraciones siguientes. Vemos que la convergencia se produjo antes del final del período de calentamiento y es estable a lo largo de las iteraciones porque las cuatro cadenas se mezclaron bien.

A partir de la imagen 4 (y de la ausencia de advertencia por parte de stan ), podemos concluir que el modelo convergió.

Evaluación de los parámetros estimados

La estadística bayesiana considera los parámetros como estocásticos, por lo que estima una distribución para cada parámetro. En la práctica, stan almacena las estimaciones de los parámetros para cada iteración posterior al calentamiento para cada cadena, lo que nos ofrece 500 x 4 estimaciones de la distribución a posteriori.

Podemos extraer del objeto stanfit un resumen de la distribución estimada utilizando la función summary:

# summarise estimated parameters
estimated <- summary(fit, pars=pars)$summary

estimated %>% kbl() %>% kable_minimal()
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu 7.538576 0.0429317 1.615445 4.395117 6.471525 7.538612 8.600434 10.75700 1415.884 1.001074
sigma 49.685819 0.0280643 1.131912 47.445472 48.895331 49.656773 50.447527 51.91045 1626.738 1.001501

A continuación, podemos comparar los parámetros \(\hat\mu\) y \(\hat\sigma\) con la media y la desviación típica observadas de \(\hat{y}\), así como con el valor verdadero, mediante la función stan_plot.

# plot estimated parameters
mean_pop <- mean(data$y)
sd_pop <- sd(data$y)


mu_plot <- stan_plot(fit, 'mu',fill_color='orange')+
  annotate('segment',x=mean_pop, xend=mean_pop, 
           y=0.7, yend=1.2,col='grey40', size=1)+
  annotate('text',x=mean_pop, 
           y=1.5, col='grey40',label= paste0('Observed\n',round(mean_pop,2)), fontface =2, size=4.5)+
    annotate('segment',x=5, xend=5, 
           y=0.7, yend=1.2,col='grey10', size=1)+
    annotate('text',x=5, 
           y=1.5, col='grey10',label= paste0('True\n',5), fontface =2, size=4.5)+
    annotate('text',x=estimated['mu', 'mean'], 
           y=1.5, col='orange',label= paste0('Estimated\n',round(estimated['mu', 'mean'],2)), fontface =2, size=4.5)
sigma_plot <- stan_plot(fit, 'sigma', fill_color='orange')+
  annotate('segment',x=sd_pop, xend=sd_pop, 
           y=0.7, yend=1.2,col='grey40', size=1)+
  annotate('text',x=sd_pop, 
           y=1.5, col='grey40',
           label= paste0('Observed\n', round(sd_pop,2)), fontface =2, size=4.5)+
    annotate('segment',x=50, xend=50, 
           y=0.7, yend=1.2,col='grey10', size=1)+
    annotate('text',x=50, 
           y=1.5, col='grey10',label= paste0('True\n',50), fontface =2, size=4.5)+
    annotate('text',x=estimated['sigma', 'mean'], 
           y=1.5, col='orange',label= paste0('Estimated\n',round(estimated['sigma', 'mean'],2)), fontface =2, size=4.5)

gridExtra::grid.arrange(mu_plot, sigma_plot, nrow=2)

Vemos que (1) la media y la desviación típica observadas están dentro de los intervalos de credibilidad del 95 % de los parámetros estimados, (2) la media y la desviación típica verdaderas están dentro de los intervalos de credibilidad del 95 % de los parámetros estimados,

La estructura del modelo se ajusta a los datos observados y consigue aproximarse al verdadero proceso de generación de datos.

Probemos con datos reales

Descarguemos los datos que vamos a modelizar. Pertenece al material complementario del documento principal que describe los modelos de población ascendentes de WorldPop (Leasure et al. 2020).

# download  tutorial data
download.file(
  "https://www.pnas.org/highwire/filestream/949050/field_highwire_adjunct_files/1/pnas.1913050117.sd01.xls",
  'tutorials/data/nga_demo_data.xls',
  method='libcurl',
  mode='wb'
)

Los datos

Los datos consisten en encuestas de hogares que recopilaron información sobre la población total en 1141 conglomerados en 15 de los 37 estados de Nigeria durante 2016 y 2017. Los grupos variaban ligeramente de tamaño, pero todos contaban con aproximadamente 3 hectáreas. Estos grupos de población se eligieron de forma aleatoria y sus límites se trazaron a partir de imágenes de satélite de alta resolución. Las encuestas se describen con más detalle en Leasure et al. (2020) y Weber et al. (2018). En la imagen 5 se muestra la ubicación de los sitios de las encuestas.

Microcensus surveys mapped as the number of survey locations within a 20km grid cell. Region groupings are shaded and numbered R1 - R11. Source: Leasure et al. (2020).

Imagen 5: Encuestas microcensales calculadas como el número de sitios de encuestas en una cuadrícula de 20 km. Los grupos de regiones están sombreados y enumeradas de R1 a R11. Fuente: Leasure et al. (2020).

El mapa de la imagen 5 muestra algunas características clave del diseño de la muestra aleatoria estratificada:

  • Solo se tomaron muestras de algunos Estados

  • Pero se tomaron muestras de al menos 1 Estado por “región”

  • Dentro de los Estados, los sitios se muestrearon aleatoriamente dentro de los tipos de asentamientos

Veamos los atributos en la tabla:

#load data
data <- readxl::read_excel(here('tutorials/data/nga_demo_data.xls'))
# create unique cluster id
data <- data %>% 
  mutate(
    id= paste0('cluster_',(1:n())), # compute cluster id
    pop_density = N/A # compute population density
  )
data %>% select(id,N, A) %>% head(10) %>% kbl() %>% kable_minimal()
id N A
cluster_1 547 3.036998
cluster_2 803 2.989821
cluster_3 750 3.078278
cluster_4 281 3.019307
cluster_5 730 3.025204
cluster_6 529 3.090072
cluster_7 505 3.007513
cluster_8 402 3.019307
cluster_9 388 3.202116
cluster_10 900 3.054689

Cada fila es un sitio de encuesta, con recuentos de población (N) y la superficie asentada (A) en hectáreas.

Variable de respuesta: el recuento de la población

Queremos modelizar la distribución del recuento de la población en cada sitio de encuesta:

# plot population count
ggplot(data, aes(x=N))+
  geom_histogram(bins=50)+
  theme_minimal()+
  theme(axis.text.y = element_blank())+
  labs(title = "", y='', x='Observed population count')+
  geom_vline(xintercept = mean(data$N), color='orange', size=1)+
  annotate('text', x=500, y=25, label=paste0('Observed mean: ', round(mean(data$N))), colour='orange', angle=90)
Observed population count distribution at survey sites

Imagen 6: Distribución del recuento de la población observada en los sitios de encuesta

Observa la gran variación en el recuento de población por sitio de encuesta, con un máximo de 2370 personas.

Modelización del recuento de población: una distribución de Poisson

El recuento de población es, por definición, una variable discreta y positiva. Una buena distribución es la distribución de Poisson.

\[ population \sim Poisson( \lambda ) \]

El DAG correspondiente muestra la interacción entre el conteo de población y el parámetro del modelo:

A continuación, tenemos que definir la distribución a priori de \(\lambda\), que corresponde a la media de la distribución de Poisson. Elegiremos una distribución a priori relativamente desinformada basada en nuestra comprensión de los datos.

Sabemos que el recuento medio de la población observada es de unos 450 habitantes por grupo. Configuramos la distribución a priori para que \(\lambda\) siga una variable uniforme entre 0 y 3000 para asegurar un parámetro positivo a la vez que reducimos la información dada al proceso de estimación y, por lo tanto, la tendencia introducida.

# define lambda prior
data$prior_lambda <- runif(nrow(data), min=0, max=3000)

ggplot(data)+
  geom_histogram(aes(x=N, after_stat(density)), bins=50, fill='grey30')+
  geom_density(aes(x=prior_lambda), colour='#00BFC4', size=1)+
  theme_minimal()+
  annotate('text', y=0.0005, x=2000, label=paste0('Uniform(0,3000)'), colour='#00BFC4', size=5)+
    theme(axis.text.y = element_blank())+
  labs(y='', x='Observed population count')
Lambda prior distribution

Imagen 7: Distribución a priori lambda

El modelo final es:

\[\begin{equation} population \sim Poisson( \lambda ) \\ \\ \lambda \sim Uniform( 0, 3000 )\tag{3} \end{equation}\]

Aplicación del modelo

Para estimar \(\lambda\), escribiremos nuestro primer modelo de población en stan.

// Model 1: Population count as a Poisson distribution 
data{
  int<lower=0> n; // number of microcensus clusters
  int<lower=0> population[n]; // count of people
}
parameters{
  // rate
  real<lower=0> lambda; 
}

model{
  // population totals
  population ~ poisson(lambda);
  // rate
  lambda ~ uniform(0, 3000);
}

Declaramos la variable de entrada population como un número entero porque nuestros datos de población son recuentos y la configuramos para que sea positiva. Definimos \(\lambda\) como un número real positivo.

Almacenamos este modelo en tutorial1_model1.stan.

Estimación del modelo

Preparamos los datos para stan:

# prepare data for stan
stan_data <- list(
  population = data$N,
  n = nrow(data))

Mantenemos los mismos parámetros que antes para el algoritmo de la cadena de Markov y declaramos \(\lambda\) como el parámetro a controlar

# parameters to monitor
pars <- c('lambda')

¡Y ya estamos listos para ejecutar el modelo!

# mcmc
fit <- rstan::stan(file = file.path('tutorial1_model1.stan'), 
                   data = stan_data,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)

El trazado muestra un modelo que converge hacia la media observada:

traceplot(fit, inc_warmup=T)

Observa la gran variación al principio de la estimación. Esto se debe a que la distribución a priori uniforme para \(\lambda\) declara que cualquier valor entre 0 y 3000 es igualmente posible.

Evaluación de la bondad de ajuste del modelo

Parámetros estimados

Trazamos el parámetro estimado \(\hat\lambda\) para ver cómo se compara con el recuento medio de población observado.

# plot estimated parameters
mean_pop <- mean(data$N)

stan_plot(fit, 'lambda',fill_color='orange')+
  annotate('segment',x=mean_pop, xend=mean_pop, 
           y=0.7, yend=1.2,col='grey40', size=1)+
  annotate('text',x=mean_pop, 
           y=1.5, col='grey40',label= paste0('Observed average\n',round(mean_pop,2)), fontface =2, size=4.5)

La media estimada se aproxima mucho a la media observada.

Recuento de población previsto

Para ver si el modelo es coherente con las observaciones, podemos calcular el recuento de población previsto para cada sitio de encuesta. Forma parte de la comprobación predictiva posterior, que se basa en la siguiente idea: si un modelo se ajusta bien, deberíamos poder utilizarlo para generar datos que se parezcan mucho a los datos que observamos.

Es posible hacerlo en stan como parte del proceso de estimación a través del bloque de cantidades generadas.

// Model 1bis: Population count as a normal distribution with integrated predictions
...
generated quantities{
   real population_hat[n];

   for(idx in 1:n){
     population_hat[idx] = poisson_rng( lambda );
   }
}

Definimos el parámetro population_hat de manera aleatoria (representado por el sufijo rng, que corresponde a random number generator –generador de números aleatorios–) de una distribución de Poisson con el parámetro estimado \(\hat\lambda\) para cada iteración.

Ejecutamos el modelo almacenado en tutorial1_model1bis.stan:

pars <- c('lambda', 'population_hat')

# mcmc
fit_model1 <- rstan::stan(file = file.path('tutorial1_model1bis.stan'), 
                   data = stan_data,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)

Y extraemos el recuento de población previsto.

# extract predictions
predicted_pop_model1 <- as_tibble(extract(fit_model1, 'population_hat')$population_hat)

colnames(predicted_pop_model1) <- data$id

Obtenemos una tabla con 500 predicciones * 4 cadenas para cada sitio de encuesta.

predicted_pop_model1 %>% 
  mutate(iteration= paste0('iter_', 1:(iter*chains))) %>% 
  select(iteration, 1:10) %>% head(10) %>% kbl() %>% kable_minimal()
iteration cluster_1 cluster_2 cluster_3 cluster_4 cluster_5 cluster_6 cluster_7 cluster_8 cluster_9 cluster_10
iter_1 459 432 446 444 432 430 433 444 432 452
iter_2 478 432 432 470 438 456 442 479 480 435
iter_3 438 414 451 456 445 390 439 456 444 475
iter_4 463 452 503 459 455 465 425 452 431 458
iter_5 421 462 444 440 468 407 399 475 448 438
iter_6 439 414 437 450 456 418 429 454 436 438
iter_7 464 467 489 470 445 455 406 459 429 481
iter_8 447 432 452 471 432 436 428 440 472 445
iter_9 424 427 432 447 444 482 477 472 458 401
iter_10 461 439 426 474 438 425 469 445 451 431

Obtenemos así una distribución de predicción a posteriori del recuento de población para cada sitio de encuesta. La imagen 8 muestra una distribución a posteriori para el primer sitio de encuesta.

# plot posterior prediction for one cluster
ggplot(predicted_pop_model1, aes(x=cluster_1))+
  geom_density(size=1.5, color='orange')+
  theme_minimal()+
    theme(axis.text.y = element_blank())+
  labs(title = "Population prediction for cluster 1 ", y='', x='')
Example of posterior prediction distribution for one cluster

Imagen 8: Ejemplo de distribución de predicción a posteriori para un grupo

Podemos extraer la predicción media para cada sitio de encuesta y un intervalo de credibilidad del 95 %.

# summarize predictions
comparison_df <- predicted_pop_model1 %>% 
      pivot_longer(everything(),names_to = 'id', values_to = 'predicted') %>% 
      group_by(id) %>% 
      summarise(across(everything(), list(mean=~mean(.), 
                                          upper=~quantile(., probs=0.975), 
                                          lower=~quantile(., probs=0.025))))
comparison_df %>% head() %>% kbl() %>% kable_minimal()
id predicted_mean predicted_upper predicted_lower
cluster_1 445.1085 488.000 402.000
cluster_10 445.0675 488.000 404.000
cluster_100 445.1840 487.025 404.975
cluster_1000 445.6165 484.000 403.000
cluster_1001 445.2855 487.000 404.975
cluster_1002 444.4365 486.000 403.000

Observamos que todas las predicciones son muy similares. El modelo solo incluye un parámetro –lambda– y las predicciones a nivel de sitio no tienen en cuenta ninguna característica de este. Por lo tanto, las predicciones de cada sitio se extraen exactamente de la misma distribución.

Veamos el panorama global trazando el recuento de población observado frente al previsto. Un modelo perfecto vería todos los puntos en la línea 1:1.

# add observed values
comparison_df <- comparison_df %>% 
  left_join(data %>% 
              select(id, N), by = 'id')

# plot predicted vs observed
ggplot(comparison_df) +
  geom_pointrange(aes(x=N, y=predicted_mean, ymin=predicted_lower, ymax=predicted_upper
                      ),
                   fill='grey50', color='grey70', shape=21
                  )+
  geom_abline(slope=1, intercept = 0, color='orange', size=1)+
  theme_minimal()+
  labs(title = '', x='Observed population count', y='Predicted population')
Comparison between observed and predicted population count (with the Poisson model). Orange line indicates the 1:1 line

Imagen 9: Comparación entre el recuento de población observado y el previsto (con el modelo de Poisson). La línea naranja indica la línea 1:1

La imagen 9 es una gran visualización del proceso de predicción. Dado que el modelo no tiene covariables (introducidas en el tutorial 3) ni estructura jerárquica (introducida en el tutorial 2), no se introduce ninguna variación subnacional. Además, los intervalos creíbles implican pocas observaciones: pocas líneas grises se cruzan con la línea naranja. Esto indica problemas con la modelización.

Podemos calcular las métricas de bondad de ajuste para completar el cuadro:

  • El patrón, la media de los valores residuales (predicción - observación)

  • La imprecisión, desviación típica del valor residual

  • La inexactitud, media de los valores residuales absolutos

  • La proporción de observaciones que se encuentran dentro del intervalo creíble previsto

# compute goodness-of-fit metrics
comparison_df %>%
  mutate(residual = predicted_mean-N,
          in_CI = ifelse(N>predicted_lower &N<predicted_upper, T, F)) %>% 
  summarise(
    `Bias`= mean(residual),
    `Imprecision` = sd(residual),
    `Inaccuracy` = mean(abs(residual)),
    `Correct credible interval (in %)` = round(sum(in_CI)/n()*100,1)
  ) %>% 
    kbl(caption = "Poisson model goodness-of-fit metrics") %>% kable_minimal()
Table 1: Poisson model goodness-of-fit metrics
Bias Imprecision Inaccuracy Correct credible interval (in %)
0.0293225 315.8384 230.9604 10

La tabla 1 confirma que el modelo está mal especificado: solo el 10 % de las observaciones se encuentra dentro de sus intervalos de credibilidad.

Esta limitación se debe a la imposibilidad de modelizar la sobredispersión en un marco de Poisson. En efecto, \(\lambda\) define tanto la media como la varianza de una variable de Poisson.

Una fuente de sobredispersión proviene del tamaño de los grupos. Observamos recuentos de población para unidades con diferentes áreas, en particular diferentes tamaños de zonas pobladas como se muestra en la imagen 10.

#plot overdispersion
ggplot(data, aes(x=A, y=N))+
  geom_point()+
  theme_minimal()+
  labs(x='Settled area in hectares', y='Population count')
Scatterplot contrasting population count with settled area for every cluster

Imagen 10: Diagrama de dispersión que contrasta el recuento de población con la zona poblada para cada grupo

Modelización del recuento de población: un modelo de Poisson lognormal

Para incorporar la sobredispersión en nuestro modelo, descomponemos el recuento de población del siguiente modo:

\[ population = pop\_density * settled\_area \]

El lado izquierdo de la ecuación es la población, una variable discreta observada, mientras que el lado derecho se compone de dos variables continuas: settled_area es la zona poblada, una variable observada, y pop_density es la densidad de población, una variable latente.

Escribir el modelo

Formalmente, podemos reescribir la ecuación (3) de la siguiente manera:

\[ population \sim Poisson( pop\_density * settled\_area) \]

Esta formulación desvela la variable latente continua positiva libre de escala, pop_density, que puede modelizarse con su propia distribución.

Optamos por una Lognormal , que es una distribución de probabilidad continua de una variable aleatoria cuyo logaritmo se distribuye normalmente. Se caracteriza por una distribución positiva sesgada a la derecha. La distribución lognormal tiene dos parámetros: \(\mu\), su parámetro de localización que define su media, y \(\sigma\), su parámetro de escala que define su desviación estándar geométrica. Este modelo nos permite reflejar la sobredispersión a través de \(\sigma\).

Probability distribution functions of log-normal distributions Source: Krishnavedala, CC0, via Wikimedia Commons

Imagen 11: Funciones de distribución de probabilidad de las distribuciones lognormales. Fuente: Krishnavedala, CC0, vía Wikimedia Commons

Aplicado a nuestra modelización de la población arroja:

Y en el modelado formal:

\[\begin{equation} population \sim Poisson( pop\_density * settled\_area) \\ pop\_density \sim Lognormal( \mu, \sigma) \end{equation}\]

Tenga en cuenta que esta ecuación es equivalente a:

\[\begin{equation} log(pop\_density) \sim Normal( \mu, \sigma) \end{equation}\]

Bajo esta forma vemos que el modelo no es lineal, sino log-lineal. Se puede reescribir como:

\[\begin{equation} pop\_density \sim exp(Normal( \mu, \sigma)) \end{equation}\]

Definición de la distribución a priori

En la distribución lognormal, hay dos parámetros, \(\mu\) que representa la media de la densidad de población en la escala logarítmica, y \(\sigma\), la desviación estándar geométrica de la densidad de población en la escala logarítmica.

Establecemos sus variables a priori de forma similar a la anterior y obtenemos de los datos que la media logarítmica observada de la densidad de población es 4,82 y la desviación geométrica estándar logarítmica observada es 0,87.

Elegimos como variable a priori para \(\mu\) una distribución normal (5,4):

# 5 Model2: Lognormal ----

# define prior for mu
data$prior_mu <- rnorm(nrow(data), mean= 5, sd=4)

ggplot(data)+
  geom_histogram(aes(x=log(pop_density), after_stat(density)), bins=100, fill='grey30')+
  geom_density(aes(x=prior_mu), colour='#00BFC4', size=1)+
  theme_minimal()+
  annotate('text', y=0.1, x=15, label=paste0('Normal(5,4)'), colour='#00BFC4', size=5)+
    geom_vline(xintercept = median(log(data$pop_density)), color='orange', size=1)+
  annotate('text', x=5.1, y=0.3, label=paste0('Observed median: ', round(median(log(data$pop_density))), 2), colour='orange', angle=90)+
    theme(axis.text.y = element_blank())+
  labs(y='', x='Observed population density (log)')+
  xlim(-8,17)
Prior distribution for mu

Imagen 12: Distribución a priori de mu

Elegimos como variable a priori para \(\sigma\) una distribución uniforme (0,4)

# define prior for sigma
data$prior_sigma <- runif(nrow(data), min = 0, max=4)+5

ggplot(data)+
  geom_histogram(aes(x=log(pop_density), after_stat(density)), bins=100, fill='grey30')+
  geom_density(aes(x=prior_sigma), colour='#00BFC4', size=1, trim=T)+
  theme_minimal()+
  annotate('text', y=0.3, x=7.7, label=paste0('Uniform(0,4)'), colour='#00BFC4', size=5)+
    theme(axis.text.y = element_blank())+
  labs(y='', x='Observed population density (log)')
Prior distribution for mu

Imagen 13: Distribución a priori de mu

El modelo resultante puede escribirse del siguiente modo:

\[\begin{equation} population \sim Poisson( pop\_density * settled\_area) \\ pop\_density \sim Lognormal( \mu, \sigma) \\ \\ \mu \sim Normal( 5, 4 ) \\ \sigma \sim Uniform( 0, 4 )\tag{4} \end{equation}\]

Aplicación del modelo

Adaptamos el código stan al cambio de modelo que afecta a todos los bloques de código:

// Model 2: Population count as a Poisson-Lognormal distribution 
data{
  int<lower=0> n; // number of microcensus clusters
  int<lower=0> population[n]; // count of people
  vector<lower=0>[n] area; // settled area
}
parameters{
  // population density
  vector<lower=0>[n] pop_density;
  // intercept
  real mu; 
  // variance
  real<lower=0> sigma; 
}
model{
  // population totals
  population ~ poisson(pop_density .* area);
  pop_density ~ lognormal( mu, sigma );
  // intercept
  mu ~ normal(5, 4);
  // variance
  sigma ~ uniform(0, 4);
}
generated quantities{
   int<lower=0> population_hat[n];
   real<lower=0> density_hat[n];

   for(idx in 1:n){
     density_hat[idx] = lognormal_rng( mu, sigma );
     population_hat[idx] = poisson_rng(density_hat[idx] * area[idx]);
   }
}

Ten en cuenta que tenemos que utilizar el formato vector para area y pop_density. Este formateo está implícito en la definición de la operación de multiplicación para la tasa de Poisson.

Almacenamos el modelo en tutorial1_model2.stan y preparamos los datos correspondientes:

# prepare data for stan
stan_data_model2 <- list(
  population = data$N,
  n = nrow(data),
  area = data$A)

Ten en cuenta que la densidad de población no es un dato de entrada del modelo, sino una variable latente no observada que modelizamos a través de la distribución lognormal.

A continuación, declaramos los parámetros a controlar (incluido density_hat) y ejecutamos el modelo.

# set parameters to monitor
pars <- c('mu','sigma', 'population_hat', 'density_hat')

# mcmc
fit_model2 <- rstan::stan(file = file.path('tutorial1_model2.stan'), 
                   data = stan_data_model2,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)

No se emiten advertencias.

Pregunta: ¿Puedes trazar el gráfico e interpretarlo?

Haz clic para ver la solución

Trazamos el diagrama:

# plot trace
traceplot(fit_model2, c('mu', 'sigma'))

Vemos que las cadenas se han mezclado bien y que las distribuciones a posteriori de los parámetros no están limitadas por su especificación a priori.

Comparamos la densidad prevista y el recuento previsto con las observaciones:

# extract posterior predictions
predicted_pop_model2 <- as_tibble(extract(fit_model2, 'population_hat')$population_hat)

predicted_dens_model2 <- as_tibble(extract(fit_model2, 'density_hat')$density_hat)
colnames(predicted_pop_model2) <- data$id
colnames(predicted_dens_model2) <- data$id

# summarise posterior predictions
comparison_df <- rbind(predicted_dens_model2 %>% 
   pivot_longer(everything(),names_to = 'id', values_to = 'predicted') %>% 
   group_by(id) %>% 
   summarise(across(everything(), list(mean=~mean(.), 
                                       upper=~quantile(., probs=0.975), 
                                       lower=~quantile(., probs=0.025)))) %>% 
   mutate(source= 'Poisson-Lognormal model',
          type= 'Population density') %>% 
   left_join(data %>% 
               select(id, pop_density)%>% 
              rename(reference=pop_density), by = 'id'),
  predicted_pop_model2 %>% 
  pivot_longer(everything(),names_to = 'id', values_to = 'predicted') %>% 
  group_by(id) %>% 
  summarise(across(everything(), list(mean=~mean(.), 
                                      upper=~quantile(., probs=0.975), 
                                      lower=~quantile(., probs=0.025)))) %>% 
  mutate(source= 'Poisson-Lognormal model',
         type='Population count') %>% 
  left_join(data %>% 
              select(id, N) %>% 
              rename(reference=N), by = 'id'))
# plot posterior predictions
ggplot(comparison_df %>% 
         mutate(type= factor(type, levels=c('Population density', 'Population count')))) +
  geom_pointrange(aes(x=reference, y=predicted_mean, ymin=predicted_lower, ymax=predicted_upper
                      ),
                   fill='grey50', color='grey70', shape=21
                  )+
  geom_abline(slope=1, intercept = 0, color='orange', size=1)+
  theme_minimal()+
  labs(title = '', x='Observations', y='Predictions')+ 
  facet_wrap(.~type, scales = 'free')

Para la variable population_density vemos el mismo patrón de estimación que en el modelo de Poisson, es decir, una distribución a posteriori de predicción media similar para cada uno de los sitios de encuesta. En cambio, la longitud de los intervalos de confianza varía entre las predicciones de los grupos debido al término de varianza de la distribución lognormal. Además, el parámetro population_count previsto se ve influido por la variable settled_area.

Utilizamos la misma métrica basada en los valores residuales para evaluar la bondad de ajuste del modelo.

# compute goodness-of-fit metrics
comparison_df %>%
  filter(type=='Population count') %>% 
  mutate(residual = predicted_mean-reference,
          in_CI = ifelse(reference>predicted_lower &reference<predicted_upper, T, F)) %>% 
  summarise(
    `Bias`= mean(residual),
    `Imprecision` = sd(residual),
    `Inaccuracy` = mean(abs(residual)),
    `Correct credible interval (in %)` = round(sum(in_CI)/n()*100,1)
  ) %>% 
    kbl(caption = "Poisson-Lognormal model goodness-of-fit metrics") %>% kable_minimal()
Table 2: Poisson-Lognormal model goodness-of-fit metrics
Bias Imprecision Inaccuracy Correct credible interval (in %)
61.51073 338.2973 265.7858 95.2

La proporción de observaciones que caen dentro de sus intervalos de credibilidad muestra un modelo con un buen comportamiento: El 95,5 % de las observaciones se encuentra en el intervalo de credibilidad del 95 %. Por término medio, el modelo sobreestima el recuento de la población en unas 60 personas, pero esta estadística varía por término medio en 337 personas.

En este modelo no se incluyen variaciones: creamos un modelo nacional de recuento de población. Este requerirá nuevas mejoras que se introducirán en los próximos tutoriales.

Y, para ello, almacenaremos este último modelo como un archivo RDS.

saveRDS(fit_model2, 'tutorial1_model2_fit.rds')

Agradecimientos

Este tutorial ha sido redactado por Edith Darin, de WorldPop, Universidad de Southampton, y Douglas Leasure, del Leverhulme Centre for Demographic Science, Universidad de Oxford, con la supervisión de Andrew Tatem, de WorldPop, Universidad de Southampton.

Este trabajo ha sido financiado por el Fondo de Población de las Naciones Unidas (UNFPA, por sus siglas en inglés).

Licencia

El presente documento puede redistribuirse libremente bajo los términos de una licencia Creative Commons Atribución-SinDerivadas 4.0 Internacional (CC BY-ND 4.0).

Referencias

Hoffman, Matthew D y Andrew Gelman. s.f. “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo”, 31.
Leasure, Douglas R., Warren C. Jochem, Eric M. Weber, Vincent Seaman y Andrew J. Tatem. 2020. “National population mapping from sparse survey data: A Hierarchical Bayesian Modeling Framework to Account for Uncertainty”. Proceedings of the National Academy of Sciences 117 (39): 24173–79. https://doi.org/10.1073/pnas.1913050117.
Weber, Eric M., Vincent Y. Seaman, Robert N. Stewart, Tomas J. Bird, Andrew J. Tatem, Jacob J. McKee, Budhendra L. Bhaduri, Jessica J. Moehl y Andrew E. Reith. 2018. “Census-Independent Population Mapping in Northern Nigeria”. Remote Sensing of Environment 204 (January): 786-98. https://doi.org/10.1016/j.rse.2017.09.024.

  1. Más detalles aquí.↩︎