logo

# 1 Set-up ----

# load libraries
library(tidyverse) # managing data
library(ggdag) # drawing DAG
library(kableExtra) # visualising table
library(here) # handling path
library(rstan) # running Bayesian models
library(plotly) # interactive plot

# stan setup
options(mc.cores = parallel::detectCores()-1)
rstan::rstan_options(auto_write = TRUE) # speed up running time of compiled model

knitr::opts_chunk$set(fig.align = "center")

Introducción

En el tutorial 1, modelizamos el recuento de población con un modelo Poisson lognormal. Esta modelización no incluye las variaciones del recuento de población a nivel de sitio.

El tutorial 2 explorará los patrones espaciales de los recuentos de población y cómo se agrupan en el espacio.

Analizaremos la agrupación de nuestras observaciones y desentrañaremos su estructura jerárquica. A continuación, veremos cómo los modelos bayesianos ofrecen diversas herramientas de modelización para representar datos jerárquicos. Por último, probaremos diferentes estructuras de modelización jerárquica para ver cuál se ajusta mejor a nuestros datos.

Los modelos jerárquicos definen un espacio de parámetros complejo. Así pues, revisaremos la configuración previa para tener en cuenta la complejidad.

Objetivos

Comprender el concepto de observaciones agrupadas
Estimar la densidad de población en un marco sin agrupación
Estimar la densidad de población en un marco con agrupación parcial
Integrar un nivel jerárquico adicional a la modelización

Lecturas complementarias

An Introduction to Hierarchical Modeling Michael Freeman. Great visualisation of grouped data and goals of general hierarchical modelling
Bayes rules! on hierarchical model Alicia A. Johnson, Miles Ott, Mine Dogucu. Good explanation of hierarchical modelling in stan with many examples
Exchangeability and hierarchical models in Bayesian Data Analysis (p. 104), Andrew Gelman, John Carlin, Hal Stern, Donald Rubin, David Dunson y Aki Vehtari

Paquetes adicionales

Utilizaremos el paquete bayesplot que ofrece más flexibilidad en la evaluación de estimaciones bayesianas y plotly para el trazado interactivo.

install.packages('bayesplot') # additional evaluations of stan models
install.packages('plotly') # interactive plot

Encontrarás viñetas aquí: bayesplot y plotly.

Estructura jerárquica de los datos

Desentrañemos los patrones de agrupamiento en los recuentos de población. La imagen 1 muestra la variación espacial de la densidad de población en Nigeria. Elegimos representar la densidad de población (personas/hectáreas pobladas) porque es más comparable entre los distintos sitios de encuesta.

Podemos observar patrones agrupados, como altas densidades de población en el norte y en el sureste.

Imagen 1: Mapa de las densidades de población de los sitios de encuesta superpuesto a los límites de la región de Nigeria

Tracemos la distribución de la densidad de población por regiones:

# 2 Hierarchy in the data ----

library(RColorBrewer)

# prepare data
data <- readxl::read_excel(here('tutorials/data/nga_demo_data.xls'))
data <- data %>% 
  mutate(
    id = as.character(1:nrow(data)),
    pop_density=N/A
  )

# plot population density per region
ggplot(data %>% 
         group_by(
           region
         ) %>% 
         mutate(mean_popDens = mean(pop_density)) %>% 
         ungroup(), aes(fill=mean_popDens, x=pop_density, y=as.factor(region)))+
  geom_boxplot()+
  theme_minimal()+
  scale_fill_stepsn( colours = brewer.pal(6, "YlOrRd"))+
  labs(fill='Mean \npopulation \ndensity', x='Population density', y='Region')

Imagen 2: Boxplot de densidades de población por región

Además de las 11 regiones, los datos proporcionan otras dos divisiones administrativas nigerianas: estatal (15) y local (es decir, áreas de Gobierno local, 223), lo que ofrece una flexibilidad adicional para agrupar los lugares sitios de encuesta. El supuesto subyacente: la estructura administrativa de un país proporciona información sobre la variación de la densidad de población.

Otro agrupamiento clave es el atributo de type indicado en los datos. Se refiere a un mapa de asentamientos basado en imágenes de satélite (Pleiades, Airbus y WorldView2, DigitalGlobe) clasificado en ocho tipos de asentamiento representados en la imagen 3 Weber et al. (2018).

Ejemplos de los tipos residencial urbano (A-F), residencial rural (M) y no residencial (Z) para los estados de Kano y Kaduna, Nigeria; Weber et al. (2018)

Imagen 3: Ejemplos de los tipos residencial urbano (A-F), residencial rural (M) y no residencial (Z) para los estados de Kano y Kaduna, Nigeria; Weber et al. (2018)

# plot population density per settlement type
ggplot(data %>% 
         group_by(
           type
         ) %>% 
         mutate(mean_popDens = mean(pop_density)) %>% 
         ungroup(), aes(fill=mean_popDens, x=pop_density, y=as.factor(type)))+
  geom_boxplot()+
  theme_minimal()+
  scale_fill_stepsn( colours = brewer.pal(6, "YlOrRd"))+
  labs(fill='Population density \n(mean)', x='', y='Settlement type')

Imagen 4: Diagrama de caja de densidades de población por tipo de asentamiento

Vemos que el tipo de asentamiento estratifica claramente la densidad de población; por ejemplo, el tipo de asentamiento 1 tiene una densidad de población considerablemente mayor que el tipo de asentamiento 4.

Obsérvese que de los ocho tipos de asentamientos que aparecen en la imagen 3, solo cinco están realmente presentes en los datos:

En la zona de tipo no residencial (Z) no se encuestó ningún sitio.
Los tipos C, D y E se han fusionado para reforzar su prevalencia en la región del sur

Imagen completa del agrupamiento de datos

La imagen 5 muestra la estructura completa de los sitios de encuesta, con la siguiente jerarquización: tipo de asentamiento, región, Estado y área de Gobierno local. Al pasar el ratón sobre el esquema, puedes ver cuántos sitios de encuesta hay en cada grupo. Puedes acceder a la estructura detallada haciendo clic en un grupo específico.

library(plotly)
# Visualise grouping of the data
# create unique id for the nested admin level
data1 <- data %>% 
  mutate(state= paste0(state,region),
         local = paste0(state, local))

# create data for sunburst plot
d1 <- rbind(
  # first layer
  data1 %>% 
    group_by(type) %>% 
    summarise(n=n()) %>% 
    mutate(
      ids = paste0('settlement', type),
      labels = paste0('settlement <br>', type),
      parents = '') %>% 
    ungroup() %>% 
    select(ids,labels, parents,n),
  # second layer
  data1 %>% 
    group_by(type, region) %>% 
    summarise(n=n()) %>% 
    mutate(
      ids = paste('settlement', type, '-', 'region', region),
      labels = paste0('region ', region),
      parents = paste0('settlement', type))%>% 
    ungroup() %>% 
    select(ids,labels, parents,n),
  # third layer
  data1 %>% 
    group_by(type, region, state) %>% 
    summarise(n=n()) %>% 
    mutate(
      ids = paste('settlement', type, '-', 'region', region, '-', 'state', state, '-', 'region', region),
      labels = paste0('state ', state),
      parents = paste('settlement', type, '-', 'region', region))%>% 
    ungroup() %>% 
    select(ids,labels, parents,n),
  # fourth layer
  data1 %>% 
    group_by(type, region, state, local) %>% 
    summarise(n=n()) %>% 
    mutate(
      ids = paste('settlement', type, '-', 'region', region, '-', 'state', state,  '-', 'local', local),
      labels = paste0('local ', local),
      parents = paste('settlement', type, '-', 'region', region, '-', 'state', state, '-', 'region', region))%>% 
    ungroup() %>% 
    select(ids,labels, parents,n)
) %>% 
  mutate(
    hover= paste(ids, '\n sample size', n)
  )

plot_ly(d1, ids = ~ids, labels = ~labels, parents = ~parents, type = 'sunburst', 
        hovertext=~hover, insidetextorientation='radial')

Imagen 5: Estructura de agrupación completa de la encuesta (asentamiento, región, estado, local)

Ten en cuenta lo siguiente:

debido a las limitaciones del muestreo, no todos los Estados y municipios están presentes en los datos
no todas las combinaciones de grupos están representadas en los datos. Algunas regiones/Estados/áreas de Gobierno local no disponen de sitios de encuesta pertenecientes a todos los tipos de asentamiento.

La imagen 6 muestra las regiones en las que falta el tipo de asentamiento. Este análisis no nos permite determinar si esos tipos de asentamiento faltan porque no se han incluido en nuestra encuesta o porque no existen (por ejemplo, una región remota que no contiene ningún tipo de asentamiento altamente urbano).

makeSunburst2layer <- function(data){
  layers <- rbind(
  # first layer
  data %>% 
    group_by(type) %>% 
    summarise(n=sum(!is.na(N))) %>% 
    mutate(
      ids = paste0('settlement', type),
      labels = paste0('settlement <br>', type),
      parents = '') %>% 
    ungroup() %>% 
    select(ids,labels, parents,n),
  # second layer
  data %>% 
    group_by(type, region) %>% 
    summarise(n=sum(!is.na(N))) %>% 
    mutate(
      ids = paste('settlement', type, '-', 'region', region),
      labels = paste0('region ', region),
      parents = paste0('settlement', type))%>% 
    ungroup() %>% 
    select(ids,labels, parents,n)) %>%
    mutate(
      hover= paste(ids, '\n sample size', n),
      color= ifelse(n==0, 'yellow','')
    )
  
 return(layers)
}

# create missing combinations
data1_complete <- data1 %>% 
  complete(region, nesting(type)) 


plot_ly() %>% 
   add_trace(data=makeSunburst2layer(data1), 
             ids = ~ids, labels = ~labels, parents = ~parents, 
             type = 'sunburst', 
             hovertext=~hover, marker= list(colors=~color),  
             insidetextorientation='radial',
             domain = list(column = 0)) %>% 
     add_trace(data=makeSunburst2layer(data1_complete), 
             ids = ~ids, labels = ~labels, parents = ~parents, 
             type = 'sunburst', 
             hovertext=~hover, marker= list(colors=~color),  
             insidetextorientation='radial',
             domain = list(column = 1))  %>%  
  layout(
      grid = list(columns =2, rows = 1),
      margin = list(l = 0, r = 0, b = 0, t = 0))

Imagen 6: Estructura de agrupamiento de la encuesta (asentamiento, región) con las combinaciones que faltan en amarillo

Estructura jerárquica en el modelo

Cuando todas las observaciones se tratan juntas e indistintamente en el modelo, se denomina agrupación completa. Así es como modelizamos el recuento de población en el tutorial ¹. Supone una capacidad de intercambio perfecta entre las observaciones, es decir, la indexación de los grupos no influye en la probabilidad de la densidad de población. En otras palabras, supone que no se pueden ordenar ni agrupar los datos. Por lo general, cuanto menos sepamos sobre un problema, con más confianza podremos afirmar la capacidad de intercambio. La consecuencia de la agrupación completa es que se estiman parámetros únicos para todo el conjunto de datos, lo que podría producir conclusiones erróneas que son válidas a nivel global pero inválidas a nivel de grupo, al debilitar los patrones específicos de grupo (véase la falacia ecológica).

Cuando un modelo se ajusta para cada grupo de datos de forma independiente, se denomina sin agrupación y crea un modelo por grupo. Reduce drásticamente el tamaño de la muestra e imposibilita la extrapolación para un grupo que no está presente en la muestra de observación.

El objetivo de la modelización jerárquica es lograr una agrupación parcial en la que se comparta la información entre grupos, teniendo en cuenta, al mismo tiempo, la estructura jerárquica de los datos. Ofrece información sobre:

la variabilidad entre grupos: las densidades de población difieren de una región a otra
la variabilidad dentro del grupo: algunas regiones presentan una mayor variabilidad que otras (véase la región 9 del gráfico 2)

Relación entre modelos aleatorios y modelos jerárquicos: Los modelos aleatorios son cualquier modelo que asuma que los parámetros no son fijos, sino que varían entre grupos, es decir, que tienen un efecto aleatorio. Los modelos jerárquicos son modelos aleatorios que proporcionan una estructura general para extraer información de otros grupos.

Recordemos nuestro modelo anterior:

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

\(\mu\) y \(\sigma\) tienen ambos un efecto fijo: se estiman para todo el conjunto de datos a la vez.

Una intercepción jerárquica por tipo de asentamiento

Queremos diferenciar los parámetros en función de la estructura jerárquica de los datos.

Un primer parámetro es \(\mu\), es decir, la media en la escala logarítmica de la distribución de la densidad de población. Modelizar \(\mu\) jerárquicamente significa reconocer que la media difiere en función de la agrupación, por ejemplo, es diferente para el tipo de asentamiento 1 que para el tipo de asentamiento 4. A partir de ahora descompondremos \(\mu\) en un término \(\alpha\)

Veamos las distintas situaciones de estimación de \(\alpha\) por tipo de asentamiento (solo mostraremos para tres tipos de asentamiento por simplicidad):

d1 <- dagify(
    Y ~ alpha,
    outcome = 'Y'
  ) %>%
    tidy_dagitty(seed=1) %>% 
    mutate(color=c('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=3,  parse=T) +
    scale_shape_manual(values=c(15,19))+
    theme_dag()+ guides(color='none', shape='none')+ labs(title = 'Complete pooling', color='', shape='')

d2 <- dagify(
    Y ~ alpha_1+alpha_2+ alpha_3,
    outcome = 'Y'
  ) %>%
    tidy_dagitty(seed=1) %>% 
    mutate(color=c('parameter','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=3,  parse=T) +
    scale_shape_manual(values=c(15,19))+
    theme_dag()+ guides(color='none', shape='none')+ labs(title = 'No pooling', color='', shape='')

d3 <- dagify(
    Y ~ alpha_1+alpha_2+ alpha_3,
    alpha_1 ~ alpha,
    alpha_2 ~ alpha,
    alpha_3 ~ alpha,
    outcome = 'Y',
    latent = 'alpha'
  ) %>%
    tidy_dagitty(seed=7) %>% 
    mutate(color=c('parameter','parameter','parameter','parameter','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=3,  parse=T) +
    scale_shape_manual(values=c(15,19))+
    theme_dag()+ guides(color='none', shape='none')+ labs(title = 'Partial pooling', color='', shape='')

gridExtra::grid.arrange(d1,d2,d3, ncol=3)

Veremos la diferencia entre los tres modelos: agrupación completa (un \(\alpha\)), sin agrupación (3 \(\alpha\)) y agrupación parcial (3 \(\alpha\) + 1 \(\alpha\) global).

El modelo 2 del tutorial 1 se basa en la agrupación completa de los datos o efecto fijo: un \(\alpha\) único para todo tipo de asentamiento.

Sin agrupación

Implementemos un marco sin agrupación, es decir, \(\alpha\) se estimará de manera independiente para cada tipo de asentamiento \(t\). La ecuación (1) representa su modelo formal, en el cual \(\alpha\) está indexado por \(t\), el tipo de asentamiento, y cada \(\alpha_t\) tiene una distribución normal a priori idéntica pero independiente Normal (5, 4).

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

¿Cómo lo escribes en Stan? Solo mostraremos el código afectado.

// Model 1: Independent alpha by settlement type
data{
  ...
  int<lower=0> type[n]; // settlement type
  int<lower=0> ntype; // number of settlement types
}
parameters{
  ...
  // independent intercept by settlement type
  vector[ntype] alpha_t; 
}
model{
  // population totals
  ...
  pop_density ~ lognormal( alpha_t[type], sigma );
  
  // independent intercept by settlement type
  alpha_t ~ normal(5, 4);
  ...
}
generated quantities{
  ...
   for(idx in 1:n){
     density_hat[idx] = lognormal_rng( alpha_t[type[idx]], sigma );
     population_hat[idx] = poisson_rng(density_hat[idx] * area[idx]);
   }
}

En el bloque de datos, declaramos la estructura de observación del tipo de asentamiento.

En el bloque de parámetros, α se convierte en un vector de tamaño 5, el número del tipo de asentamiento en los datos.

En el bloque del modelo, elegimos la misma distribución a priori para cada αt, cinco distribuciones normales centradas en 5 y con desviación típica 4.

Observa cómo funciona la indexación en el bloque de cantidades generadas.

Una vez escrito el modelo, preparamos los datos correspondientes, agregando el tipo de asentamiento:

# 3. Hierarchical model by settlement type ----

# No-pooling
# prepare data for stan
stan_data_model1 <- list(
  population = data$N,
  n = nrow(data),
  area = data$A,
  type = data$type,
  ntype= n_distinct(data$type))

Y ejecutamos el modelo de forma similar al tutorial 1:

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

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

# mcmc
fit1 <- rstan::stan(file = file.path('tutorial2_model1.stan'), 
                   data = stan_data_model1,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)

stan no emite ninguna advertencia; el modelo ha convergido.

La imagen 7 muestra que los \(\hat\alpha_t\) t estimadas varían según el tipo de asentamiento. También vemos cómo el \(\alpha\) estimado previamente se promediaba entre los diferentes patrones.

# plot estimated parameter
stan_plot(fit1, pars='alpha_t', fill_color='orange')+
  # add alpha from tutorial 1
  geom_vline(xintercept=4.755109, size=1.5, linetype=2)+
  annotate('text', x=5, y=5.7, label="alpha estimated \nin tutorial 1")

Imagen 7: Intercepción estimada independientemente por tipo de asentamiento

Dibujamos el diagrama observado frente al previsto para ver el impacto de la estimación de \(\alpha\) por tipo de asentamiento:

# write function to extract posterior predictions and summarise them in a table
getPopPredictions <- function(model_fit, 
                              estimate='population_hat',
                              obs='N', reference_data=data){
  # extract predictions
  predicted_pop <- as_tibble(extract(model_fit, estimate)[[estimate]])
  colnames(predicted_pop) <- reference_data$id
  
  # summarise predictions
  predicted_pop <- predicted_pop %>% 
    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)))) %>% 
    left_join(reference_data %>% 
                rename('reference'=all_of(obs)) %>% 
                select(id, reference), by = 'id')%>% 
    mutate(
      residual= predicted_mean - reference,
      ci_size = predicted_upper- predicted_lower,
      estimate = estimate
      )

return(predicted_pop)
}

# plot posterior predictions
ggplot(
  rbind(
    getPopPredictions(fit1) %>% 
      mutate(type='Population count'),
    getPopPredictions(fit1, estimate = 'density_hat', obs='pop_density') %>% 
      mutate(type='Population density')
  )  %>% 
         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')

Imagen 8: Comparación de las observaciones con las predicciones del modelo 1

Vemos que la densidad de población prevista tiene ahora cinco modos, que corresponden a los cinco tipos de asentamiento.

Agrupación parcial

Probemos ahora la agrupación parcial. La agrupación parcial implica una configuración jerárquica de las distribuciones a priori, es decir, cada distribución a priori \(\alpha_t\) dependerá de una distribución a priori \(\alpha\) general. En otras palabras, \(\alpha\) define una distribución por países que obliga a cada \(\alpha_t\) a situarse en un rango similar.

\[\begin{equation} population \sim Poisson( pop\_density * settled\_area) \\ pop\_density \sim Lognormal( \alpha_t, \sigma) \\[5pt] \alpha_t \sim Normal( \alpha, \nu_\alpha ) \\[5pt] \alpha \sim Normal( 5, 10 ) \\ \nu_\alpha \sim Uniform( 0, 15 ) \\ \sigma \sim Uniform( 0, 10 )\tag{2} \end{equation}\]

Ten en cuenta lo siguiente:

en \(\alpha_t\), los efectos aleatorios se extraen cada uno de una distribución normal con los mismos hiperparámetros generales –parámetros de distribuciones a priori– \(\alpha\) y \(\nu_\alpha\).
Relajamos la restricción a priori sobre \(\alpha\) y \(\sigma\) de una distribución normal (5, 4) y uniforme (0, 4) a una normal (5, 10) y uniforme (0, 10). Los modelos jerárquicos son difíciles de estimar, ya que algunos parámetros se definen como correlacionados. Por lo tanto, ofrecemos al algoritmo más espacio para explorar.
Establecemos una distribución uniforme (0, 15) como hiperpriori –distribución a priori de un hiperparámetro– para \(\nu_\alpha\), la desviación estándar del efecto aleatorio a priori, porque las desviaciones estándar son positivas. Establecemos un rango más amplio para \(\nu_\alpha\) que para \(\sigma\) porque, por experiencia, las desviaciones típicas de los parámetros definidos jerárquicamente son más difíciles de estimar que otras desviaciones típicas (véase el punto 2).

En stan escribiremos el modo de la siguiente manera:

// Model 2: Hierarchical alpha by settlement type
parameters{
  ...
  // hierarchical intercept by settlement
  vector[ntype] alpha_t; 
  real alpha;
  real<lower=0> nu_alpha;
}
model{
  ...
  // hierarchical intercept by settlement
  alpha_t ~ normal(alpha, nu_alpha);
  alpha ~ normal(5, 10);
  nu_alpha ~ uniform(0, 15);
}

Añadimos los nuevos parámetros a los parámetros a controlar y ejecutamos el modelo:

# Partial pooling

pars <- c('alpha_t','alpha', 'nu_alpha', 'sigma', 'population_hat', 'density_hat')

# mcmc
fit2 <- rstan::stan(file = file.path('tutorial2_model2.stan'), 
                   data = stan_data_model1,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)

Ahora podemos comparar los parámetros estimados entre un marco sin agrupación y uno con agrupación parcial

# extract parameters from both models
fit2_alpha_t <- summary(fit2, pars=c('alpha_t'))$summary
fit1_alpha_t <-  summary(fit1, pars=c('alpha_t'))$summary

# plot parameters
data_plot <- rbind(
  as_tibble(fit1_alpha_t, rownames='param') %>% 
    mutate(model='No pooling'),
    as_tibble(fit2_alpha_t, rownames='param') %>% 
    mutate(model='Partial pooling')
) %>% 
  mutate(
    model = factor(model, levels=c('No pooling','Partial pooling')),
    param_ = paste(param,model),
    labels = ifelse(model=='Partial pooling', param, '')) %>% 
arrange(param_)

ggplot(data_plot, aes(mean,param_, color=model, fill=model,label=labels))+
  geom_point()+
  geom_linerange(aes(xmin=`2.5%`, xmax=`97.5%`))+
  theme_minimal()+
  labs(y='')+
    scale_y_discrete(labels=data_plot$labels)

Imagen 9: Comparación entre los valores alfa_t estimados de forma independiente y jerárquica

No hay diferencias notables entre los dos modelos, pero el marco sin agrupación tiende a sobreestimar la diferencia entre cada efecto aleatorio, mientras que el marco con agrupación parcial tiende a acercar los efectos aleatorios entre sí. La tabla 1 muestra cómo la incorporación del agrupamiento por tipo de asentamiento tiene un efecto directo en la predicción.

# Complete pooling 

# load previous model for complete pooling
fit_tuto1_model2 <- readRDS('../tutorial1/tutorial1_model2_fit.rds')

# build comprehensive dataframe
comparison_df <- rbind(
 getPopPredictions(fit1) %>% 
   mutate(model='No pooling'),
 getPopPredictions(fit2) %>% 
   mutate(model='Partial pooling'),
 getPopPredictions(fit_tuto1_model2) %>% 
   mutate(model='Complete pooling'))


# compute goodness-of-fit metrics
comparison_df %>% group_by(model) %>% 
  summarise( `Bias`= mean(residual),
    `Inaccuracy` = mean(abs(residual)),
        `Imprecision` = sd(residual)
) %>%  kbl(caption = 'Goodness-of-fit metrics comparison between complete pooling, no pooling, and partial pooling') %>% kable_minimal()

Table 1: Goodness-of-fit metrics comparison between complete pooling, no pooling, and partial pooling
model	Bias	Inaccuracy	Imprecision
Complete pooling	61.22336	265.6971	337.8884
No pooling	45.99138	234.5601	307.6078
Partial pooling	45.16813	234.2322	307.2914

La estimación de \(\alpha\) por tipo de asentamiento mejora el modelo, como demuestran los valores residuales más bajos (patrón e imprecisión) y la menor variación (inexactitud). El modelo jerárquico sí reduce ligeramente el patrón, pero no de forma llamativa, lo que indica que el número de sitios de encuesta por tipo de asentamiento era lo suficientemente grande como para estimar el \(\alpha_t\) de forma independiente, en un marco sin agrupación.

Sin embargo, en la próxima sección veremos la segunda ventaja de la modelización jerárquica.

Una intercepción jerárquica por tipo de asentamiento y región

Si se tiene en cuenta la estratificación por tipo de asentamiento del recuento de población, se obtienen predicciones más precisas y se reduce la incertidumbre de las predicciones (véase la tabla 1).

Podemos incluir más grupos en la modelización, como divisiones administrativas.

Incluyamos la región en la modelización:

\[\begin{equation} population \sim Poisson( pop\_density * settled\_area) \\ pop\_density \sim Lognormal( \alpha_{t,r}, \sigma) \\[10pt] \alpha_{t,r} \sim Normal(\alpha_t, \nu_t) \\ \alpha_t \sim Normal( \alpha, \nu) \\[10pt] \alpha \sim Normal( 5, 10 ) \\ \nu \sim Uniform( 0, 15 ) \\ \nu_t \sim Uniform( 0, 15 ) \\ \sigma \sim Uniform( 0, 10 ) \end{equation}\]

Vemos que las distribuciones a priori para \(\alpha_{t,r}\), la media por región y por tipo de asentamiento, dependen de dos parámetros globales, \(\alpha_t\) y \(\nu_t\), que se estiman por tipo de asentamiento. Las distribuciones a priori para \(\alpha_t\), la media por tipo de asentamiento, dependen de dos parámetros globales, \(\alpha\) y \(\nu\), que se estiman a nivel nacional.

Aplicación del modelo

El modelo se traduce en stan de la siguiente manera:

// Model 3: Hierarchical alpha by settlement type and region

data{
  ...
  int<lower=1> nregion; //number of regions
  int<lower=1,upper=nregion> region[n]; // region
}
parameters{
  ...
  // hierarchical intercept by settlement and region
  real alpha;
  
  vector[ntype] alpha_t; 
  vector[nregion] alpha_t_r[ntype];
  
  real<lower=0> nu_alpha;
  real<lower=0> nu_alpha_t;
}
transformed parameters{
  vector[n] pop_density_median;
  
  for(idx in 1:n){
    pop_density_median[idx] = alpha_t_r[type[idx],region[idx]];
  }
}
model{
  pop_density ~ lognormal(pop_density_median, sigma );

  // hierarchical intercept by settlement and region
  alpha ~ normal(5, 10);
  nu_alpha ~ uniform(0, 15);
  nu_alpha_t ~ uniform(0, 15);
  
  alpha_t ~ normal(alpha, nu_alpha);
  
  for(t in 1:ntype){
    alpha_t_r[t,] ~ normal(alpha_t[t], nu_alpha_t); 
  }
  ...
}
generated quantities{
  ...
  density_hat[idx] = lognormal_rng( alpha_t_r[type[idx], region[idx]], sigma );
}

Observa el nuevo bloque de parámetros transformados: es útil mantener ordenado el bloque de tu modelo con solo las relaciones estocásticas.

Añadimos la indexación de la región a la lista de datos preparada para stan:

# 4. Hierarchical model by settlement type and region ----

# prepare data for stan
stan_data_model3 <- list(
  population = data$N,
  n = nrow(data),
  area = data$A,
  type = data$type,
  ntype= n_distinct(data$type),
  region = data$region,
  nregion = n_distinct(data$region)
  )

Y ejecuta el modelo con los nuevos parámetros a controlar:

pars <- c('alpha','alpha_t','alpha_t_r','nu_alpha', 'sigma','population_hat',  'density_hat')

# mcmc
fit3 <- rstan::stan(file = file.path('tutorial2_model3.stan'), 
                   data = stan_data_model3,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)

## Warning in validityMethod(object): The following variables have undefined
## values: population_hat[1],The following variables have undefined
## values: population_hat[2],The following variables have undefined
## values: population_hat[3],The following variables have undefined
## values: population_hat[4],The following variables have undefined
## values: population_hat[5],The following variables have undefined
## values: population_hat[6],The following variables have undefined
## values: population_hat[7],The following variables have undefined
## values: population_hat[8],The following variables have undefined
## values: population_hat[9],The following variables have undefined values:
## population_hat[10],The following variables have undefined values:
## population_hat[11],The following variables have undefined values:
## population_hat[12],The following variables have undefined values:
## population_hat[13],The following variables have undefined values:
## population_hat[14],The following variables have undefined values:
## population_hat[15],The following variables have undefined values:
## population_hat[16],The following variables have undefined values:
## population_hat[17],The following variables have undefined values:
## population_hat[18],The following variables have undefined values:
## population_hat[19],The following variables have undefined values:
## population_hat[20],The following variables have undefined values:
## population_hat[21],The following variables have undefined values:
## population_hat[22],The following variables have undefined values:
## population_hat[23],The following variables have undefined values:
## population_hat[24],The following variables have undefined values:
## population_hat[25],The following variables have undefined values:
## population_hat[26],The following variables have undefined values:
## population_hat[27],The following variables have undefined values:
## population_hat[28],The following variables have undefined values:
## population_hat[29],The following variables have undefined values:
## population_hat[30],The following variables have undefined values:
## population_hat[31],The following variables have undefined values:
## population_hat[32],The following variables have undefined values:
## population_hat[33],The following variables have undefined values:
## population_hat[34],The following variables have undefined values:
## population_hat[35],The following variables have undefined values:
## population_hat[36],The following variables have undefined values:
## population_hat[37],The following variables have undefined values:
## population_hat[38],The following variables have undefined values:
## population_hat[39],The following variables have undefined values:
## population_hat[40],The following variables have undefined values:
## population_hat[41],The following variables have undefined values:
## population_hat[42],The following variables have undefined values:
## population_hat[43],The following variables have undefined values:
## population_hat[44],The following variables have undefined values:
## population_hat[45],The following variables have undefined values:
## population_hat[46],The following variables have undefined values:
## population_hat[47],The following variables have undefined values:
## population_hat[48],The following variables have undefined values:
## population_hat[49],The following variables have undefined values:
## population_hat[50],The following variables have undefined values:
## population_hat[51],The following variables have undefined values:
## population_hat[52],The following variables have undefined values:
## population_hat[53],The following variables have undefined values:
## population_hat[54],The following variables have undefined values:
## population_hat[55],The following variables have undefined values:
## population_hat[56],The following variables have undefined values:
## population_hat[57],The following variables have undefined values:
## population_hat[58],The following variables have undefined values:
## population_hat[59],The following variables have undefined values:
## population_hat[60],The following variables have undefined values:
## population_hat[61],The following variables have undefined values:
## population_hat[62],The following variables have undefined values:
## population_hat[63],The following variables have undefined values:
## population_hat[64],The following variables have undefined values:
## population_hat[65],The following variables have undefined values:
## population_hat[66],The following variables have undefined values:
## population_hat[67],The following variables have undefined values:
## population_hat[68],The following variables have undefined values:
## population_hat[69],The following variables have undefined values:
## population_hat[70],The following variables have undefined values:
## population_hat[71],The following variables have undefined values:
## population_hat[72],The following variables have undefined values:
## population_hat[73],The following variables have undefined values:
## population_hat[74],The following variables have undefined values:
## population_hat[75],The following variables have undefined values:
## population_hat[76],The following variables have undefined values:
## population_hat[77],The following variables have undefined values:
## population_hat[78],The following variables have undefined values:
## population_hat[79],The following variables have undefined values:
## population_hat[80],The following variables have undefined values:
## population_hat[81],The following variables have undefined values:
## population_hat[82],The following variables have undefined values:
## population_hat[83],The following variables have undefined values:
## population_hat[84],The following variables have undefined values:
## population_hat[85],The following variables have undefined values:
## population_hat[86],The following variables have undefined values:
## population_hat[87],The following variables have undefined values:
## population_hat[88],The following variables have undefined values:
## population_hat[89],The following variables have undefined values:
## population_hat[90],The following variables have undefined values:
## population_hat[91],The following variables have undefined values:
## population_hat[92],The following variables have undefined values:
## population_hat[93],The following variables have undefined values:
## population_hat[94],The following variables have undefined values:
## population_hat[95],The following variables have undefined values:
## population_hat[96],The following variables have undefined values:
## population_hat[97],The following variables have undefined values:
## population_hat[98],The following variables have undefined values:
## population_hat[99],The following variables have undefined values:
## population_hat[100],The following variables have undefined values:
## population_hat[101],The following variables have undefined values:
## population_hat[102],The following variables have undefined values:
## population_hat[103],The following variables have undefined values:
## population_hat[104],The following variables have undefined values:
## population_hat[105],The following variables have undefined values:
## population_hat[106],The following variables have undefined values:
## population_hat[107],The following variables have undefined values:
## population_hat[108],The following variables have undefined values:
## population_hat[109],The following variables have undefined values:
## population_hat[110],The following variables have undefined values:
## population_hat[111],The following variables have undefined values:
## population_hat[112],The following variables have undefined values:
## population_hat[113],The following variables have undefined values:
## population_hat[114],The following variables have undefined values:
## population_hat[115],The following variables have undefined values:
## population_hat[116],The following variables have undefined values:
## population_hat[117],The following variables have undefined values:
## population_hat[118],The following variables have undefined values:
## population_hat[119],The following variables have undefined values:
## population_hat[120],The following variables have undefined values:
## population_hat[121],The following variables have undefined values:
## population_hat[122],Th

Stan nos dice que “las siguientes variables tienen valores indefinidos: population_hat[1]”. Esto se debe a un desbordamiento de enteros: el límite del valor entero es 2^31 - 1. Para asegurarnos de no sobrepasar este límite, stan restringe el valor de la tasa de Poisson a 1,07374e+09.

Sin embargo, este problema solo se produce en la predicción y durante el período de calentamiento. Cuando las cadenas alcanzan la convergencia, la tasa de Poisson pasa a ser razonable y el parámetro population_hat queda definido. Aquí presentamos el ejemplo del proceso de predicción del grupo 2:

traceplot(fit3, 'density_hat[2]', inc_warmup=T)

Imagen 10: Traceplot para el clúster 2

Parámetros estimados

Veamos la \(\hat\alpha_{t,r}\) estimada. Dado que hay 1 \(\hat\alpha\) , 5 \(\hat\alpha\) t y 5x11 \(\hat\alpha_{t,r}\), nos centraremos solo en dos tipos de asentamientos en concreto, el 1 y el 4. Para ello, recuerda cómo se construyen los nombres de los parámetros:

alpha_t es un vector de dimensión 5, de modo que alpha_t[1] representa \(\hat{\alpha_1}\), la estimación asociada al asentamiento 1
alpha_t_r es un vector de dimensión 5*11, de modo que alpha_t[1,1] representa \(\hat{\alpha}_{1,1}\), la estimación asociada al asentamiento 1 y a la región 1.

alpha_4_r <- paste0('alpha_t_r[4,', 1:11, ']')
alpha_1_r <- paste0('alpha_t_r[1,', 1:11, ']')

stan_plot(fit3, pars=c('alpha','alpha_t[1]','alpha_t[4]', alpha_1_r, alpha_4_r), 
          fill_color='orange')

Imagen 11: Alfa estimado por región para los asentamientos 1 y 4

Vemos que:

\(\hat\alpha_1\) tiene una distribución estimada que se superpone a los dos tipos de asentamiento y un intervalo de confianza más amplio para dar cuenta de toda la incertidumbre
\(\hat\alpha_1\) 1 y \(\hat\alpha_1\) 2 presentan dos patrones distintivos
los dos patrones a nivel de asentamiento ocultan disparidades regionales

Una cuarta observación: analiza detenidamente las regiones 4 y 10, es decir, \(\hat\alpha_{1,4}\), \(\hat\alpha_{1,10}\), \(\hat\alpha_{2,4}\) y \(\hat\alpha_{2,10}\). Vemos que los intervalos de confianza son mayores que para cualquier otra región.

Esto se debe a que esas dos regiones solo comprenden un tipo de asentamiento (5) en los datos (véase la imagen 6) y, por lo tanto, ni el tipo de asentamiento 1 ni el tipo de asentamiento 4. Por lo tanto, no hay datos para informar a esos \(\alpha_{t,r}\), pero, debido a la jerarquía, se pueden estimar a partir del parámetro global \(\alpha_t\). Y, en efecto, vemos que sus medias estimadas coinciden estrechamente con los respectivos parámetros \(\hat\alpha_t\).

Esta característica de la estructura jerárquica es una gran ventaja en comparación con la estructura sin agrupación. Al predecir el recuento de población en toda Nigeria, podríamos encontrar un tipo de asentamiento en una región que no estuviera representado en nuestro conjunto de datos, normalmente el tipo 1, 2, 3 o 4 en la región 4 o 10. En una estructura sin agrupación no tendremos ningún parámetro estimado porque los modelos se estiman independientemente por combinación de niveles para una región de tipo x. En una estructura de agrupación parcial o un modelo jerárquico, estamos cubiertos por el \(\alpha_t\) global. Obviamente, habrá más incertidumbre en la predicción para esta zona específica porque no se han podido utilizar datos para ajustar el modelo.

Esta capacidad de un modelo jerárquico se mantiene también si no se ha muestreado un grupo entero, normalmente si no se ha incluido una región en la encuesta. Sin embargo, significa que esta región que falta se estimará con observaciones de las otras regiones muestreadas. Si la región que falta tiene características muy diferentes, el modelo no podrá captarlas.

Por último, podemos ver el efecto de una jerarquía de dos niveles en las predicciones del modelo para cada sitio de encuesta.

# build comprehensive dataframe
comparison_df <- rbind(
  comparison_df,
 getPopPredictions(fit3) %>% 
   mutate(model='Hierarchical - settlement, region'))

# compute goodness-of-fit metrics
comparison_df %>% 
  mutate(model =ifelse(model=='Partial pooling', 'Hierarchical - settlement', model)) %>% 
  filter(grepl('Hierarchical', model))  %>% group_by(model) %>% 
  summarise( `Bias`= mean(residual),
    `Inaccuracy` = mean(abs(residual)),
        `Imprecision` = sd(residual)
) %>%  kbl(caption = 'Goodness-of-metrics comparison of the hierarchical models') %>% kable_minimal()

Table 2: Goodness-of-metrics comparison of the hierarchical models
model	Bias	Inaccuracy	Imprecision
Hierarchical - settlement	45.16813	234.2322	307.2914
Hierarchical - settlement, region	35.19189	200.0586	284.0779

Guardaremos el modelo 3 para utilizarlo como comparación en el tutorial 2.

# save model
saveRDS(fit3, 'tutorial2_model3_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

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.

Un excelente ejemplo del concepto de intercambiabilidad se ofrece en el libro de análisis bayesiano de datos, propuesto como lectura complementaria.↩︎

Modelización de la población para apoyar los censos

Tutorial 2: ¿Cómo modelizar las variaciones espaciales a gran escala?

Elaborado el 24/1/2022