Clicky

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

Introducción

Del tutorial 1 al tutorial 3, creamos un modelo de población ascendente, modelizando ladrillo a ladrillo. En el tutorial 1, modelizamos el recuento de población como un compuesto Poisson lognormal, teniendo en cuenta la zona poblada. En el tutorial 2, añadimos una intercepción aleatoria jerárquica por tipo de asentamiento y región, que diferencia la estimación de parámetros por la agrupación natural de las observaciones. En el tutorial 3, modelizamos las variaciones a pequeña escala de la densidad de población añadiendo una regresión lineal basada en covariables geoespaciales. De este modo, cubrimos todos los bloques de modelización fundamentales del modelo de población ascendente WorldPop, descrito en Leasure et al. (2020).

En este tutorial abordaremos el diagnóstico avanzado de modelos para una evaluación completa de la bondad de ajuste. Una vez comprobado que el modelo supera satisfactoriamente los distintos diagnósticos, procederemos a predecir el recuento de población de cada cuadrícula de la zona de estudio. La predicción de la población será la oportunidad para hablar de la incertidumbre de las estimaciones.

Objetivos

  1. Analizar los problemas de convergencia
  2. Comprender la comprobación predictiva a posteriori
  3. Comprender la validación cruzada y el sobreajuste
  4. Predecir la población de cada cuadrícula del sitio de estudio
  5. Visualizar la incertidumbre de las predicciones
  6. Predicción agregada en unidades espaciales personalizadas

Lecturas complementarias

Paquetes adicionales

Utilizaremos algunos paquetes adicionales, principalmente para la manipulación de SIG:

  • raster para predecir la cuadrícula de población en formato ráster

  • sf para agregar la predicción en unidad espacial personalizada

  • tmap para trazar datos espaciales

  • RColorBrewer para utilizar una bonita paleta de colores en los diagramas

install.packages('raster')
install.packages('RColorBrewer')
install.packages('sf')
install.packages('tmap')

Encontrarás viñetas aquí: raster, RColorBrewer, sf and tmap.

Diagnóstico avanzado de modelos

Antes de predecir el recuento de población de toda la zona de estudio, queremos asegurarnos de que el modelo es correcto. Así pues, realizamos comprobaciones adicionales: sobre la convergencia del modelo con trazados y advertencias de stan, sobre los parámetros estimados con una comprobación predictiva a posteriori y sobre el recuento de población previsto en el sitio de estudio con un enfoque de validación cruzada.

Trabajaremos con el primer modelo de población relevante que creamos en el tutorial 1, es decir, el modelo “básico” de Poisson lognormal:

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

Haremos algunas configuraciones que ahora deberían resultarnos familiares.

# 1. Model diagnostics ----

# load 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
  )
# mcmc settings
chains <- 4
iter <- 500
seed <- 1789

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

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

Evaluación de la convergencia

En esta sección exploraremos los modelos que tienen problemas de convergencia y descubriremos cómo evaluar de dónde provienen y cómo solucionarlos.

Veamos primero el impacto de la longitud de calentamiento, la fase inicial de las simulaciones en la que las secuencias se acercan a la masa de la distribución.

Elegimos un calentamiento de 20 iteraciones (en lugar de 250).

# 1.1 Short warmup period ----

# set short warmup
warmup <- 20

# mcmc
fit_lowWarmup <- rstan::stan(file = file.path('../tutorial1/tutorial1_model2.stan'), 
                          data = stan_data,
                          iter = warmup + iter, 
                          chains = chains,
                          warmup = warmup, 
                          pars = pars,
                          seed = seed)
## Warning: There were 1097 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems

Stan emite varias advertencias. Puedes explorar su significado aquí.

En este caso, el trazado es muy informativo:

traceplot(fit_lowWarmup, c('mu', 'sigma'), inc_warmup=T)

Los trazados indican que se sigue explorando después de la zona sombreada, es decir, una vez finalizado el periodo de calentamiento.

A continuación, examinaremos un problema de convergencia más grave. Supongamos que elegimos valores a priori inadecuados para nuestro modelo. Normalmente podemos declarar una distribución uniforme con límites más estrictos para \(\sigma\):

\[\begin{equation} \sigma 〜 Uniform( 0, 1 ) \end{equation}\]

Almacenamos el modelo en tutorial4/tutorial1_model2_wrong.stan y lo estimamos.

# 1.2 Wrong prior ----

warmup <- 250

# mcmc
fit_wrong <- rstan::stan(file = file.path('tutorial1_model2_wrong.stan'), 
                   data = stan_data,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)
## Warning: There were 1981 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 4.4, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess

Muchas métricas detectan un problema de convergencia. La más preocupante son las transiciones divergentes que suelen estar relacionadas con problemas en la propia estructura del modelo y no solo dentro del proceso de estimación.

traceplot(fit_wrong, c('mu', 'sigma'))

El trazado de \(\mu\) indica que las cadenas no se mezclaron. Además, el trazado de \(\sigma\) muestra claramente el problema de limitar la estimación con un tope de 1 que la restringe en exceso.

Comprobación de la predicción a posteriori

Volvamos a un modelo correctamente ajustado.

# 2.3 Predicted posterior check ----

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

Como se ha mostrado en el párrafo anterior, la estimación bayesiana depende en gran medida de la elección de los valores a priori correctos para los parámetros. Esos valores a priori deben reflejar el conocimiento experto sobre la cantidad de interés que se confronta, a continuación, con las observaciones. El modelo a priori define el espacio para la estimación a posteriori. Por lo tanto, debemos asegurarnos de que el apoyo a priori no restringe excesivamente el apoyo a posteriori. Esto se denomina comprobación predictiva a posteriori.

Primero tenemos que extraer la distribución a posteriori estimada del objeto stanfit. Lo hacemos utilizando la función extract de rstan.

# extract estimated mu
mus <- data.frame(
    # extract parameter
    posterior=extract(fit, pars='mu')$mu,
    parameter = 'mu'
    )
glimpse(mus)
## Rows: 2,000
## Columns: 2
## $ posterior <dbl> 4.765670, 4.735604, 4.788244, 4.726230, 4.746974, 4.783040, ~
## $ parameter <chr> "mu", "mu", "mu", "mu", "mu", "mu", "mu", "mu", "mu", "mu", ~

mus contiene la mu estimada para cada iteración posterior al calentamiento (500) de cada cadena de Markov (4).

Queremos recrear la distribución a priori. En el modelo 2 del tutorial 3, utilizamos la siguiente distribución a priori: \[\mu \sim Normal(5,4)\] a partir de esta distribución a priori, el mismo número de valores que tenemos para la distribución a posteriori (es decir, 2000).

# retrieve stan parameter
post_warmup_draws <- iter - warmup

# simulate from the prior distribution
mus$prior <- rnorm(chains * post_warmup_draws, 5, 4)

Para \(\sigma\), creamos un marco de datos similar a mus que contrasta su distribución a priori con su distribución a posteriori y lo unimos a mus.

# build dataframe with parameters distribution comparison
parameters <-  rbind(
  # alpha
  mus ,
  #sigma
  data.frame(
      posterior=extract(fit, pars='sigma')$sigma,
      prior = runif(chains * post_warmup_draws, 0, 4),
      parameter = 'sigma')
  )  %>% 
  pivot_longer(-parameter,
    names_to = 'distribution',
    values_to = 'value'
  )

# plot distribution comparison for both parameter
ggplotly(ggplot(parameters, aes(x=value, color=distribution, fill=distribution))+
  geom_density()+
  theme_minimal()+
    facet_wrap(.~parameter, scales = 'free'))

Imagen 1: Comprobaciones predictivas a posteriori de alfa y sigma

La imagen 1 muestra una comprobación predictiva a posteriori tanto para \(\mu\) como para \(\sigma\). Compara la distribución a priori que fijamos en el modelo y la distribución a posteriori estimada. No dudes en ampliar el zoom, ya que el rango de valores es muy diferente la distribución a priori y la distribución a posteriori.

Confirmamos que la distribución a priori no restringe falsamente la distribución a posteriori, con un rango muy amplio de valores posibles. La única restricción que podemos observar es el límite inferior de \(\sigma\) en 0, pero esto es normal porque una varianza tiene que ser positiva.

Validación cruzada de las predicciones del modelo

Tras comprobar la estructura del modelo, volvemos a examinar los parámetros de bondad de ajuste basados en las predicciones de población.

En los tutoriales 1 a 3, comprobamos la bondad de ajuste de las predicciones modelizadas en la muestra, es decir, prediciendo la población y comparándola con las observaciones de los lugares de estudio que ya se utilizaron para ajustar el modelo.

Las comprobaciones dentro de la muestra no tienen en cuenta el sobreajuste del modelo. El sobreajuste se produce cuando el modelo no solo reproduce la información de las observaciones, sino también el ruido. Queremos que el modelo sea capaz de predecir con exactitud la población también para nuevos lugares de estudio desconocidos. Este escenario del mundo real puede reproducirse conservando un porcentaje de las observaciones para el ajuste del modelo y dejando el resto de los lugares de estudio para la predicción, un mecanismo denominado validación cruzada.

Primero tenemos que dividir nuestros datos en dos: un conjunto de entrenamiento (70 %), que se utilizará para ajustar el modelo y un conjunto de predicción (30 %), que se utilizará para la predicción del modelo:

# 2.4 Cross validation of predictions ----

set.seed(2004)
# sample observations
train_size <-  round(nrow(data)*0.7)
train_idx <- sample(1:nrow(data), size = train_size, replace=F)

# build train datasets
train_data <- data[train_idx,]

# build test datasets
test_data <- data[-train_idx,]

En stan, la validación cruzada puede realizarse simultáneamente con el ajuste del modelo proporcionando el conjunto de datos de prueba al bloque de modelización de generated quantities. Como hay que definir todas las variables para cada conjunto, se alarga sustancialmente el código.

data{
  
  int<lower=0> n_train; // number of microcensus clusters in training set
  int<lower=0> n_test; // number of microcensus clusters in predicting set

  int<lower=0> population_train[n_train]; // count of people
  
  vector<lower=0>[n_train] area_train; // settled area in training
  vector<lower=0>[n_test] area_test; // settled area in testing
}

parameters{
  // population density
  vector<lower=0>[n_train] pop_density_train;
  
  //intercept
  real mu;

  // variance
  real<lower=0> sigma; 
}

model{
  
  // population totals
  population_train ~ poisson(pop_density_train .* area_train);
  pop_density_train ~ lognormal( mu, sigma );
  
  //  intercept
  mu ~ normal(5, 4);
  
  // variance
  sigma ~ uniform(0, 4);
}

generated quantities{
  
  int<lower=-0> population_hat[n_test];
  real<lower=0> density_hat[n_test];

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

Preparamos los datos para stan y ejecutamos el modelo:

# prepare data
stan_data_xval <- list(
  population_train = train_data$N,
  n_train = nrow(train_data),
  n_test = nrow(test_data),

  area_train = train_data$A,
  area_test = test_data$A,

  seed=seed
)

# mcmc setting
pars <- c('mu','sigma', 'population_hat',  'density_hat')

# mcmc
fit_xval <- rstan::stan(file = file.path('tutorial1_model2xval.stan'), 
                   data = stan_data_xval,
                   iter = warmup + iter, 
                   chains = chains,
                   warmup = warmup, 
                   pars = pars,
                   seed = seed)

Ahora tenemos acceso a dos conjuntos de predicciones: uno, fit, de la predicción dentro de la muestra y otro, fit_xval, de la predicción de validación cruzada.

Comparamos la bondad de las métricas de los dos modelos:

# predict population
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)))
      ) %>% 
    # merge with observations
    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)
}

# build comparison dataframe between in-sample and xvalidation
comparison_df <- rbind(
 getPopPredictions(fit) %>% 
   mutate(Model='In-sample'),
 getPopPredictions(fit_xval, reference_data = test_data) %>% 
   mutate(Model='Cross-validation'))

# 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-metrics computed in-sample vs cross-validation') %>% kable_minimal()
Table 1: Goodness-of-metrics computed in-sample vs cross-validation
Model Bias Inaccuracy Imprecision
Cross-validation 96.40268 263.2451 313.2041
In-sample 61.22336 265.6971 337.8884

La tabla 1 muestra que todas las métricas muestran un peor ajuste: el patrón, la imprecisión y la inexactitud aumentan. Sin embargo, se mantienen en el mismo rango, lo que parece indicar que nuestro modelo no está sobreajustado.

Una buena práctica consiste en evaluar siempre el ajuste del modelo en el conjunto de pruebas para evitar producir métricas de bondad de ajuste demasiado optimistas.

Predicción de la población en cuadrículas

En esta sección trataremos las predicciones del modelo para toda la zona de estudio a nivel de celda de cuadrícula. La población puede predecirse para cualquier unidad espacial siempre que 1. todas las covariables y los datos de entrada puedan producirse para la unidad espacial 2. la unidad espacial de predicción no sea demasiado diferente en tamaño de la unidad espacial de ajuste, de forma que la relación estimada entre las covariables y la densidad de población pueda seguir aplicándose.

Predeciremos la población para una cuadrícula de 100 m x 100 m y utilizaremos nuestro mejor modelo de población: el modelo 2 del tutorial 3.

# 3. Population prediction ----

tutorial3_model2_fit <- readRDS("../tutorial3/tutorial3_model2_fit.rds")

Esto requiere extraer las distribuciones de los parámetros estimados: \(\hat\alpha_{t,r}, \hat\sigma, \hat\beta^{fixed}, \hat\beta^{random}_t\). A continuación, estos parámetros estimados se combinan con los datos ráster de entrada siguiendo este modelo modificado para la predicción:

\[\begin{equation} population\_predicted 〜 Poisson( \hat{pop\_density} * settled\_area) \\ \hat{pop\_density} 〜 Lognormal(\hat\mu, \: \hat\sigma) \\ \hat\mu = \hat\alpha_{t,r} + \hat\beta^{random}_t X^{random} + \hat\beta^{fixed} X^{fixed} \tag{2} \end{equation}\]

\(X\) y \(settled\_area\) son los datos de entrada a nivel de celda de cuadrícula:

La predicción del modelo basada en celdas de cuadrícula requiere, por lo tanto, 10 cuadrículas de entrada:

  • zona poblada

  • tipo de asentamiento

  • región administrativa

  • las seis covariables

  • una cuadrícula maestra que define las celdas de la cuadrícula para las que predecir la población. Estas celdas cuadriculadas (1) se encuentran dentro de la zona de estudio, (2) se consideran pobladas

Lamentablemente, no se nos permite redistribuir el mapa de asentamientos utilizado en Nigeria, versión 1.2, que fue proporcionado por el Oak Ridge National Laboratory (2018). Por lo tanto, proporcionaremos el código para mostrar todos los pasos de la predicción de población, pero no podrá ejecutarse.

Nos centraremos en la región 8, que solo contiene 671 110 celdas, en lugar de las 166 412 498 celdas de todo el país.

Preparar los datos para la predicción

Para predecir la población en cuadrículas, convertimos el conjunto de datos ráster en una tabla: en las filas, las celdas de la cuadrícula y en las columnas, los valores ráster. No necesitamos conservar la estructura espacial del ráster para la predicción del modelo. Cada celda de la cuadrícula puede estimarse de forma independiente siempre que conozcamos su región administrativa, el tipo de asentamiento y los atributos de las covariables.

library(raster)
wd_path <- 'Z:/Projects/WP517763_GRID3/Working/NGA/ed/bottom-up-tutorial/'

# function that loads the raster and transforms it into an array
getRasterValues <- function(x){
  return( raster(paste0(wd_path,x))[])
}

# get a list of the rasters
rasterlist <- list.files(wd_path,
                         pattern='.tif', full.names = F)
# apply function to raster list
raster_df <- as_tibble(sapply(rasterlist, getRasterValues))

La tabla ráster contiene 671 110 filas y 10 columnas. el 95 % de las filas son NA porque no se consideran pobladas.

settled <- raster_df %>% filter(mastergrid==1)
settled %>% head() %>% kbl() %>% kable_minimal()
mastergrid settled_area x1 x2 x3 x4 x5 x6 gridcell_id settlementType
1 0.3764444 -0.0819475 0.0126183 4.663749 -0.1662194 -0.1876594 -0.1262890 47 5
1 0.0235278 -0.0629407 0.0189274 4.664463 -0.1640482 -0.1852137 -0.1265523 48 5
1 0.0411736 -0.0036624 0.0000000 4.661014 -0.1278646 -0.1549067 -0.1385387 90 5
1 0.0176458 0.0634254 0.0000000 4.662555 -0.1018127 -0.1279219 -0.1385805 91 5
1 0.0117639 0.0478863 1.2744479 4.871751 0.9074216 1.0802944 0.2204351 208 5
1 0.0999931 0.0948045 1.3028392 4.877525 0.8835372 1.0512134 0.2195443 209 5

Necesitamos escalar las covariables con los factores de escala calculados durante la etapa de ajuste.

#load scaling factor
scale_factor <- read_csv('tutorials/data/scale_factor.csv')

scaled_covs <- settled %>% 
  # subset covs column
  dplyr::select(starts_with('x'), gridcell_id) %>% 
  
  # convert table to long format
  pivot_longer(-gridcell_id, names_to='cov', values_to = 'value') %>%
  
  # add scaling factor
  left_join(scale_factor) %>% 
  
  # scale covs
  mutate(value_scaled= (value-cov_mean)/cov_sd) %>% 
  dplyr::select(gridcell_id, cov, value_scaled) %>% 
  
  # convert table back to wide format
  pivot_wider(names_from = cov, values_from = value_scaled)

# replace the covs with their scaled version
raster_df <- raster_df %>% dplyr::select(-starts_with('x')) %>% 
  left_join(scaled_covs)

Extracción de los parámetros estimados

El segundo paso consiste en extraer las distribuciones de parámetros estimados. Utilizamos el modelo ajustado con el conjunto de datos de las observaciones completas y no el de la validación cruzada para obtener la máxima información.

#extract betas
beta_fixed <- t(rstan::extract(model_fit, pars='beta_fixed')$beta_fixed)
beta_random <- t(rstan::extract(model_fit, pars='beta_random')$beta_random)
#extract alphas
alphas <- rstan::extract(model_fit, pars='alpha_t_r')$alpha_t_r
#extract sigmas
sigmas <-  rstan::extract(model_fit, pars='sigma')$sigma

Primero nos centramos en producir \(\hat\beta^{fixed} X^{fixed}\) a partir de la ecuación (2), las covariables modelizadas con un efecto fijo.

\(\hat\beta^{fixed}\) es una matriz de 2000 x 5 para las 5 covariables y 2000 estimaciones. Lo transpusimos para multiplicarlo por \(X^{fixed}\), una matriz de 310 512 X 5 de covariables para cada celda asentada.

# extract covariates modelled with a fixed effect
cov_fixed <- settled %>% 
  dplyr::select(x1:x3, x5:x6) %>% 
  as.matrix()

cov_fixed <- cov_fixed %*% beta_fixed

A continuación, producimos \(\hat\beta^{random}_t X^{random}\), las covariables modelizadas con un efecto aleatorio.

\(\hat\beta^{random}_t\) es una matriz de 2000 x 5 para los 5 tipos de asentamiento. Tenemos que asociar cada celda de la cuadrícula con el parámetro \(\hat\beta^{random}_t\) correcto basado en el tipo de asentamiento de la celda de la cuadrícula y, después, multiplicar por \(X^{random}\) los valores de las covariables.

beta_random <- as_tibble(beta_random)
# add settlement type to beta random
beta_random$settlementType <- 1:5

# extract covariates modelled with a random effect
cov_random <- settled %>% 
  # subset random covariate and settlement type
  dplyr::select(settlementType,x4) %>% 
  # associate correct estimated beta_t
  left_join(beta_random) %>% 
  # multiply cov by slope
  mutate(
    across(starts_with('V'), ~ .x*x4)
  ) %>% 
  # keep only the estimates
  dplyr::select(-settlementType, -x4) %>% 
  as.matrix()

Finalmente, debemos asociar el \(\hat\alpha_{t,8}\) correcto a cada celda de la cuadrícula.

\(\hat\alpha_{t,8}\) tiene un formato similar a \(\hat\beta^{random}_t\). Procederemos de la misma manera.

# subset alpha_t for region 8
alpha_t_8 <- as_tibble(t(alphas[,,8]))
# assign settlement type
alpha_t_8$settlementType <- 1:5

alpha_predicted <- settled %>% 
  dplyr::select(settlementType) %>% 
  left_join(alpha_t_8) %>% 
  dplyr::select(-settlementType) %>% 
  as.matrix()

Finalmente, podemos calcular \(\hat\mu\) como la suma de \(\hat\alpha_{t,r}\), \(\hat\beta^{random}_t X^{random}\) y \(\hat\beta^{fixed} X^{fixed}\)

# sum mu components
mu_predicted <- alpha_predicted + cov_fixed + cov_random

Predicción del recuento de población en cada cuadrícula

Para estimar la población en celdas de la cuadrícula, necesitamos simular a partir de la ecuación (2) con los \(\hat\mu\) y \(\hat\sigma\) estimados.

predictions <- as_tibble(mu_predicted) %>% 
  # add grid cell id and the settled area to mu
  mutate(
      gridcell_id= settled$gridcell_id,
      settled_area = settled$settled_area
    )  %>% 
  # long format
  pivot_longer(
    -c(gridcell_id, settled_area), 
    names_to = 'iterations', values_to = 'mu_predicted') %>%
  mutate(
    # add sigma iterations for every grid cell
    sigma=rep(sigmas, nrow(mu_predicted)),
    # draw density for a log normal
    density_predicted = rlnorm(n(), mu_predicted, sigma),
    # draw population count from a Poisson
    population_predicted = rpois(n(), density_predicted*settled_area)
    ) %>% 
  dplyr::select(-mu_predicted,-sigma, -density_predicted) %>% 
  # convert it back to grid cell x iterations matrix
  pivot_wider(-c(iteration, population_predicted),
      names_from = iteration,
      values_from = population_predicted
    )

predictions contiene 2000 predicciones de población para cada celda de la cuadrícula poblada. Veamos qué aspecto tiene para las 5 celdas de cuadrícula pobladas:

predictions[1:5, ] %>% dplyr::select(gridcell_id, everything()) %>% 
  kbl() %>% kable_minimal()  %>% scroll_box(width = "100%")
gridcell_id V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20 V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34 V35 V36 V37 V38 V39 V40 V41 V42 V43 V44 V45 V46 V47 V48 V49 V50 V51 V52 V53 V54 V55 V56 V57 V58 V59 V60 V61 V62 V63 V64 V65 V66 V67 V68 V69 V70 V71 V72 V73 V74 V75 V76 V77 V78 V79 V80 V81 V82 V83 V84 V85 V86 V87 V88 V89 V90 V91 V92 V93 V94 V95 V96 V97 V98 V99 V100 V101 V102 V103 V104 V105 V106 V107 V108 V109 V110 V111 V112 V113 V114 V115 V116 V117 V118 V119 V120 V121 V122 V123 V124 V125 V126 V127 V128 V129 V130 V131 V132 V133 V134 V135 V136 V137 V138 V139 V140 V141 V142 V143 V144 V145 V146 V147 V148 V149 V150 V151 V152 V153 V154 V155 V156 V157 V158 V159 V160 V161 V162 V163 V164 V165 V166 V167 V168 V169 V170 V171 V172 V173 V174 V175 V176 V177 V178 V179 V180 V181 V182 V183 V184 V185 V186 V187 V188 V189 V190 V191 V192 V193 V194 V195 V196 V197 V198 V199 V200 V201 V202 V203 V204 V205 V206 V207 V208 V209 V210 V211 V212 V213 V214 V215 V216 V217 V218 V219 V220 V221 V222 V223 V224 V225 V226 V227 V228 V229 V230 V231 V232 V233 V234 V235 V236 V237 V238 V239 V240 V241 V242 V243 V244 V245 V246 V247 V248 V249 V250 V251 V252 V253 V254 V255 V256 V257 V258 V259 V260 V261 V262 V263 V264 V265 V266 V267 V268 V269 V270 V271 V272 V273 V274 V275 V276 V277 V278 V279 V280 V281 V282 V283 V284 V285 V286 V287 V288 V289 V290 V291 V292 V293 V294 V295 V296 V297 V298 V299 V300 V301 V302 V303 V304 V305 V306 V307 V308 V309 V310 V311 V312 V313 V314 V315 V316 V317 V318 V319 V320 V321 V322 V323 V324 V325 V326 V327 V328 V329 V330 V331 V332 V333 V334 V335 V336 V337 V338 V339 V340 V341 V342 V343 V344 V345 V346 V347 V348 V349 V350 V351 V352 V353 V354 V355 V356 V357 V358 V359 V360 V361 V362 V363 V364 V365 V366 V367 V368 V369 V370 V371 V372 V373 V374 V375 V376 V377 V378 V379 V380 V381 V382 V383 V384 V385 V386 V387 V388 V389 V390 V391 V392 V393 V394 V395 V396 V397 V398 V399 V400 V401 V402 V403 V404 V405 V406 V407 V408 V409 V410 V411 V412 V413 V414 V415 V416 V417 V418 V419 V420 V421 V422 V423 V424 V425 V426 V427 V428 V429 V430 V431 V432 V433 V434 V435 V436 V437 V438 V439 V440 V441 V442 V443 V444 V445 V446 V447 V448 V449 V450 V451 V452 V453 V454 V455 V456 V457 V458 V459 V460 V461 V462 V463 V464 V465 V466 V467 V468 V469 V470 V471 V472 V473 V474 V475 V476 V477 V478 V479 V480 V481 V482 V483 V484 V485 V486 V487 V488 V489 V490 V491 V492 V493 V494 V495 V496 V497 V498 V499 V500 V501 V502 V503 V504 V505 V506 V507 V508 V509 V510 V511 V512 V513 V514 V515 V516 V517 V518 V519 V520 V521 V522 V523 V524 V525 V526 V527 V528 V529 V530 V531 V532 V533 V534 V535 V536 V537 V538 V539 V540 V541 V542 V543 V544 V545 V546 V547 V548 V549 V550 V551 V552 V553 V554 V555 V556 V557 V558 V559 V560 V561 V562 V563 V564 V565 V566 V567 V568 V569 V570 V571 V572 V573 V574 V575 V576 V577 V578 V579 V580 V581 V582 V583 V584 V585 V586 V587 V588 V589 V590 V591 V592 V593 V594 V595 V596 V597 V598 V599 V600 V601 V602 V603 V604 V605 V606 V607 V608 V609 V610 V611 V612 V613 V614 V615 V616 V617 V618 V619 V620 V621 V622 V623 V624 V625 V626 V627 V628 V629 V630 V631 V632 V633 V634 V635 V636 V637 V638 V639 V640 V641 V642 V643 V644 V645 V646 V647 V648 V649 V650 V651 V652 V653 V654 V655 V656 V657 V658 V659 V660 V661 V662 V663 V664 V665 V666 V667 V668 V669 V670 V671 V672 V673 V674 V675 V676 V677 V678 V679 V680 V681 V682 V683 V684 V685 V686 V687 V688 V689 V690 V691 V692 V693 V694 V695 V696 V697 V698 V699 V700 V701 V702 V703 V704 V705 V706 V707 V708 V709 V710 V711 V712 V713 V714 V715 V716 V717 V718 V719 V720 V721 V722 V723 V724 V725 V726 V727 V728 V729 V730 V731 V732 V733 V734 V735 V736 V737 V738 V739 V740 V741 V742 V743 V744 V745 V746 V747 V748 V749 V750 V751 V752 V753 V754 V755 V756 V757 V758 V759 V760 V761 V762 V763 V764 V765 V766 V767 V768 V769 V770 V771 V772 V773 V774 V775 V776 V777 V778 V779 V780 V781 V782 V783 V784 V785 V786 V787 V788 V789 V790 V791 V792 V793 V794 V795 V796 V797 V798 V799 V800 V801 V802 V803 V804 V805 V806 V807 V808 V809 V810 V811 V812 V813 V814 V815 V816 V817 V818 V819 V820 V821 V822 V823 V824 V825 V826 V827 V828 V829 V830 V831 V832 V833 V834 V835 V836 V837 V838 V839 V840 V841 V842 V843 V844 V845 V846 V847 V848 V849 V850 V851 V852 V853 V854 V855 V856 V857 V858 V859 V860 V861 V862 V863 V864 V865 V866 V867 V868 V869 V870 V871 V872 V873 V874 V875 V876 V877 V878 V879 V880 V881 V882 V883 V884 V885 V886 V887 V888 V889 V890 V891 V892 V893 V894 V895 V896 V897 V898 V899 V900 V901 V902 V903 V904 V905 V906 V907 V908 V909 V910 V911 V912 V913 V914 V915 V916 V917 V918 V919 V920 V921 V922 V923 V924 V925 V926 V927 V928 V929 V930 V931 V932 V933 V934 V935 V936 V937 V938 V939 V940 V941 V942 V943 V944 V945 V946 V947 V948 V949 V950 V951 V952 V953 V954 V955 V956 V957 V958 V959 V960 V961 V962 V963 V964 V965 V966 V967 V968 V969 V970 V971 V972 V973 V974 V975 V976 V977 V978 V979 V980 V981 V982 V983 V984 V985 V986 V987 V988 V989 V990 V991 V992 V993 V994 V995 V996 V997 V998 V999 V1000 V1001 V1002 V1003 V1004 V1005 V1006 V1007 V1008 V1009 V1010 V1011 V1012 V1013 V1014 V1015 V1016 V1017 V1018 V1019 V1020 V1021 V1022 V1023 V1024 V1025 V1026 V1027 V1028 V1029 V1030 V1031 V1032 V1033 V1034 V1035 V1036 V1037 V1038 V1039 V1040 V1041 V1042 V1043 V1044 V1045 V1046 V1047 V1048 V1049 V1050 V1051 V1052 V1053 V1054 V1055 V1056 V1057 V1058 V1059 V1060 V1061 V1062 V1063 V1064 V1065 V1066 V1067 V1068 V1069 V1070 V1071 V1072 V1073 V1074 V1075 V1076 V1077 V1078 V1079 V1080 V1081 V1082 V1083 V1084 V1085 V1086 V1087 V1088 V1089 V1090 V1091 V1092 V1093 V1094 V1095 V1096 V1097 V1098 V1099 V1100 V1101 V1102 V1103 V1104 V1105 V1106 V1107 V1108 V1109 V1110 V1111 V1112 V1113 V1114 V1115 V1116 V1117 V1118 V1119 V1120 V1121 V1122 V1123 V1124 V1125 V1126 V1127 V1128 V1129 V1130 V1131 V1132 V1133 V1134 V1135 V1136 V1137 V1138 V1139 V1140 V1141 V1142 V1143 V1144 V1145 V1146 V1147 V1148 V1149 V1150 V1151 V1152 V1153 V1154 V1155 V1156 V1157 V1158 V1159 V1160 V1161 V1162 V1163 V1164 V1165 V1166 V1167 V1168 V1169 V1170 V1171 V1172 V1173 V1174 V1175 V1176 V1177 V1178 V1179 V1180 V1181 V1182 V1183 V1184 V1185 V1186 V1187 V1188 V1189 V1190 V1191 V1192 V1193 V1194 V1195 V1196 V1197 V1198 V1199 V1200 V1201 V1202 V1203 V1204 V1205 V1206 V1207 V1208 V1209 V1210 V1211 V1212 V1213 V1214 V1215 V1216 V1217 V1218 V1219 V1220 V1221 V1222 V1223 V1224 V1225 V1226 V1227 V1228 V1229 V1230 V1231 V1232 V1233 V1234 V1235 V1236 V1237 V1238 V1239 V1240 V1241 V1242 V1243 V1244 V1245 V1246 V1247 V1248 V1249 V1250 V1251 V1252 V1253 V1254 V1255 V1256 V1257 V1258 V1259 V1260 V1261 V1262 V1263 V1264 V1265 V1266 V1267 V1268 V1269 V1270 V1271 V1272 V1273 V1274 V1275 V1276 V1277 V1278 V1279 V1280 V1281 V1282 V1283 V1284 V1285 V1286 V1287 V1288 V1289 V1290 V1291 V1292 V1293 V1294 V1295 V1296 V1297 V1298 V1299 V1300 V1301 V1302 V1303 V1304 V1305 V1306 V1307 V1308 V1309 V1310 V1311 V1312 V1313 V1314 V1315 V1316 V1317 V1318 V1319 V1320 V1321 V1322 V1323 V1324 V1325 V1326 V1327 V1328 V1329 V1330 V1331 V1332 V1333 V1334 V1335 V1336 V1337 V1338 V1339 V1340 V1341 V1342 V1343 V1344 V1345 V1346 V1347 V1348 V1349 V1350 V1351 V1352 V1353 V1354 V1355 V1356 V1357 V1358 V1359 V1360 V1361 V1362 V1363 V1364 V1365 V1366 V1367 V1368 V1369 V1370 V1371 V1372 V1373 V1374 V1375 V1376 V1377 V1378 V1379 V1380 V1381 V1382 V1383 V1384 V1385 V1386 V1387 V1388 V1389 V1390 V1391 V1392 V1393 V1394 V1395 V1396 V1397 V1398 V1399 V1400 V1401 V1402 V1403 V1404 V1405 V1406 V1407 V1408 V1409 V1410 V1411 V1412 V1413 V1414 V1415 V1416 V1417 V1418 V1419 V1420 V1421 V1422 V1423 V1424 V1425 V1426 V1427 V1428 V1429 V1430 V1431 V1432 V1433 V1434 V1435 V1436 V1437 V1438 V1439 V1440 V1441 V1442 V1443 V1444 V1445 V1446 V1447 V1448 V1449 V1450 V1451 V1452 V1453 V1454 V1455 V1456 V1457 V1458 V1459 V1460 V1461 V1462 V1463 V1464 V1465 V1466 V1467 V1468 V1469 V1470 V1471 V1472 V1473 V1474 V1475 V1476 V1477 V1478 V1479 V1480 V1481 V1482 V1483 V1484 V1485 V1486 V1487 V1488 V1489 V1490 V1491 V1492 V1493 V1494 V1495 V1496 V1497 V1498 V1499 V1500 V1501 V1502 V1503 V1504 V1505 V1506 V1507 V1508 V1509 V1510 V1511 V1512 V1513 V1514 V1515 V1516 V1517 V1518 V1519 V1520 V1521 V1522 V1523 V1524 V1525 V1526 V1527 V1528 V1529 V1530 V1531 V1532 V1533 V1534 V1535 V1536 V1537 V1538 V1539 V1540 V1541 V1542 V1543 V1544 V1545 V1546 V1547 V1548 V1549 V1550 V1551 V1552 V1553 V1554 V1555 V1556 V1557 V1558 V1559 V1560 V1561 V1562 V1563 V1564 V1565 V1566 V1567 V1568 V1569 V1570 V1571 V1572 V1573 V1574 V1575 V1576 V1577 V1578 V1579 V1580 V1581 V1582 V1583 V1584 V1585 V1586 V1587 V1588 V1589 V1590 V1591 V1592 V1593 V1594 V1595 V1596 V1597 V1598 V1599 V1600 V1601 V1602 V1603 V1604 V1605 V1606 V1607 V1608 V1609 V1610 V1611 V1612 V1613 V1614 V1615 V1616 V1617 V1618 V1619 V1620 V1621 V1622 V1623 V1624 V1625 V1626 V1627 V1628 V1629 V1630 V1631 V1632 V1633 V1634 V1635 V1636 V1637 V1638 V1639 V1640 V1641 V1642 V1643 V1644 V1645 V1646 V1647 V1648 V1649 V1650 V1651 V1652 V1653 V1654 V1655 V1656 V1657 V1658 V1659 V1660 V1661 V1662 V1663 V1664 V1665 V1666 V1667 V1668 V1669 V1670 V1671 V1672 V1673 V1674 V1675 V1676 V1677 V1678 V1679 V1680 V1681 V1682 V1683 V1684 V1685 V1686 V1687 V1688 V1689 V1690 V1691 V1692 V1693 V1694 V1695 V1696 V1697 V1698 V1699 V1700 V1701 V1702 V1703 V1704 V1705 V1706 V1707 V1708 V1709 V1710 V1711 V1712 V1713 V1714 V1715 V1716 V1717 V1718 V1719 V1720 V1721 V1722 V1723 V1724 V1725 V1726 V1727 V1728 V1729 V1730 V1731 V1732 V1733 V1734 V1735 V1736 V1737 V1738 V1739 V1740 V1741 V1742 V1743 V1744 V1745 V1746 V1747 V1748 V1749 V1750 V1751 V1752 V1753 V1754 V1755 V1756 V1757 V1758 V1759 V1760 V1761 V1762 V1763 V1764 V1765 V1766 V1767 V1768 V1769 V1770 V1771 V1772 V1773 V1774 V1775 V1776 V1777 V1778 V1779 V1780 V1781 V1782 V1783 V1784 V1785 V1786 V1787 V1788 V1789 V1790 V1791 V1792 V1793 V1794 V1795 V1796 V1797 V1798 V1799 V1800 V1801 V1802 V1803 V1804 V1805 V1806 V1807 V1808 V1809 V1810 V1811 V1812 V1813 V1814 V1815 V1816 V1817 V1818 V1819 V1820 V1821 V1822 V1823 V1824 V1825 V1826 V1827 V1828 V1829 V1830 V1831 V1832 V1833 V1834 V1835 V1836 V1837 V1838 V1839 V1840 V1841 V1842 V1843 V1844 V1845 V1846 V1847 V1848 V1849 V1850 V1851 V1852 V1853 V1854 V1855 V1856 V1857 V1858 V1859 V1860 V1861 V1862 V1863 V1864 V1865 V1866 V1867 V1868 V1869 V1870 V1871 V1872 V1873 V1874 V1875 V1876 V1877 V1878 V1879 V1880 V1881 V1882 V1883 V1884 V1885 V1886 V1887 V1888 V1889 V1890 V1891 V1892 V1893 V1894 V1895 V1896 V1897 V1898 V1899 V1900 V1901 V1902 V1903 V1904 V1905 V1906 V1907 V1908 V1909 V1910 V1911 V1912 V1913 V1914 V1915 V1916 V1917 V1918 V1919 V1920 V1921 V1922 V1923 V1924 V1925 V1926 V1927 V1928 V1929 V1930 V1931 V1932 V1933 V1934 V1935 V1936 V1937 V1938 V1939 V1940 V1941 V1942 V1943 V1944 V1945 V1946 V1947 V1948 V1949 V1950 V1951 V1952 V1953 V1954 V1955 V1956 V1957 V1958 V1959 V1960 V1961 V1962 V1963 V1964 V1965 V1966 V1967 V1968 V1969 V1970 V1971 V1972 V1973 V1974 V1975 V1976 V1977 V1978 V1979 V1980 V1981 V1982 V1983 V1984 V1985 V1986 V1987 V1988 V1989 V1990 V1991 V1992 V1993 V1994 V1995 V1996 V1997 V1998 V1999 V2000
47 375 114 124 45 143 189 51 298 137 105 508 299 69 33 81 404 61 240 159 78 103 75 642 458 115 488 61 89 148 377 230 121 54 132 368 53 80 55 91 50 136 186 421 115 201 53 184 177 741 86 77 102 65 161 235 161 135 234 93 179 281 79 218 144 77 69 145 186 97 39 89 187 110 94 236 103 257 179 108 40 46 138 157 62 67 42 76 90 95 123 254 164 81 108 79 147 315 124 92 205 127 162 124 93 76 178 131 164 104 82 134 62 93 204 107 91 167 300 47 117 90 40 276 107 210 68 35 218 203 46 143 198 413 374 50 62 190 157 57 90 129 214 480 154 135 141 75 218 221 45 342 129 134 113 159 88 70 117 145 9 75 93 201 77 52 298 116 64 441 40 190 376 78 139 265 302 206 208 89 141 32 83 175 159 81 175 74 150 120 45 86 99 394 68 221 261 233 196 30 211 153 200 372 65 24 170 97 158 246 60 176 231 162 655 105 121 108 136 92 43 87 58 57 65 84 82 103 108 127 39 175 153 59 37 87 119 46 181 280 73 178 62 95 117 78 203 137 260 181 35 125 151 390 45 30 196 85 74 66 72 33 95 166 222 140 194 176 109 172 160 53 64 137 369 83 209 134 57 298 154 62 112 120 66 151 343 154 96 119 131 295 137 161 178 107 53 102 103 188 142 36 310 248 192 88 131 121 64 268 336 107 68 80 156 50 135 199 90 97 40 311 91 272 184 76 42 67 134 118 26 155 53 90 241 116 211 83 36 378 205 303 144 108 99 228 73 132 62 111 21 137 71 132 46 220 51 170 154 181 192 94 32 95 84 187 67 104 58 86 42 141 53 94 189 45 259 151 66 80 145 163 596 162 90 94 190 72 138 30 294 231 190 127 80 222 163 90 125 137 27 362 81 99 134 34 84 283 48 91 38 521 183 84 143 157 82 107 73 159 316 118 136 83 201 698 41 66 82 448 120 59 33 58 278 112 127 211 404 386 169 225 98 176 175 301 142 344 190 82 135 66 84 149 314 218 95 91 76 114 110 111 37 59 199 133 208 253 63 44 45 64 402 68 78 37 98 363 239 34 105 168 166 41 84 36 145 146 55 124 81 56 97 122 147 83 75 123 227 203 108 183 188 179 38 180 84 208 294 45 312 481 384 74 96 120 40 93 59 112 156 135 74 577 257 169 199 105 87 126 166 32 34 158 65 89 253 85 158 81 83 55 134 95 260 73 208 68 358 75 185 201 40 122 232 101 56 225 114 85 209 120 210 112 131 132 143 91 248 84 174 58 108 96 45 88 111 58 83 184 134 176 62 127 137 218 20 182 99 76 105 125 109 130 99 158 73 150 135 323 123 140 225 63 136 69 65 75 127 486 101 58 45 86 448 402 93 199 162 116 151 217 125 104 371 116 118 151 97 99 40 71 85 191 155 80 236 220 26 335 38 43 112 234 201 455 184 135 90 166 273 159 134 92 128 157 760 180 182 458 46 248 136 48 175 95 66 93 50 91 112 96 128 176 318 62 101 116 163 135 56 227 108 430 63 217 274 274 68 60 125 314 111 144 135 51 147 202 93 166 113 48 161 155 127 205 100 67 64 50 157 30 78 217 69 334 344 27 47 332 308 207 87 196 116 42 67 270 127 76 181 41 73 441 534 106 90 115 180 237 121 201 345 184 116 87 234 422 43 261 24 39 280 265 214 133 38 125 117 271 36 272 59 225 205 45 157 175 156 62 71 255 270 160 280 364 120 38 89 213 122 97 110 149 382 224 50 150 165 123 57 125 146 261 99 136 58 109 73 250 175 61 321 77 142 47 63 129 91 158 42 119 140 149 462 121 95 101 424 77 127 246 216 465 217 164 173 534 136 181 244 79 345 196 144 131 343 113 132 521 66 79 420 155 101 84 205 58 148 160 145 256 123 102 222 103 126 126 459 124 120 121 69 315 262 298 98 67 197 154 265 109 141 74 138 157 81 33 115 201 172 195 200 65 75 35 99 86 49 166 211 128 32 78 112 83 153 168 101 1083 48 19 118 124 92 90 27 89 31 52 54 84 30 41 311 139 123 186 109 139 628 49 210 184 50 86 197 64 115 295 200 22 40 34 257 294 152 450 35 75 288 200 178 134 303 160 126 235 117 91 83 98 78 98 136 89 118 46 331 217 72 123 66 144 335 39 258 48 67 221 70 58 33 65 152 201 101 373 141 100 237 397 89 134 199 61 91 202 88 155 120 105 152 132 83 126 298 28 159 89 33 57 122 49 478 306 172 88 92 138 134 29 99 240 107 296 223 62 41 276 89 100 75 81 88 129 68 109 144 56 282 252 93 112 124 113 112 86 15 106 68 56 75 68 183 84 81 68 249 134 209 176 132 143 71 47 44 210 85 64 61 207 77 197 79 114 66 123 31 51 139 150 74 108 154 110 149 62 157 80 260 54 53 168 298 66 228 137 320 117 143 86 129 329 364 290 130 207 104 46 82 119 46 153 151 66 119 151 109 148 93 126 87 115 182 244 393 368 78 131 273 258 117 47 198 27 108 263 459 46 55 136 88 138 149 86 420 82 162 51 95 94 213 110 93 397 258 190 67 95 114 204 108 50 185 270 77 159 85 171 107 29 126 197 323 112 56 196 103 215 56 48 82 203 119 93 123 186 131 253 44 209 112 144 96 166 84 55 185 80 117 177 273 188 122 96 123 48 81 197 231 51 93 325 180 445 250 76 78 80 164 113 157 360 258 772 202 315 520 140 181 51 154 93 130 136 112 51 339 247 43 75 240 177 201 66 34 67 143 93 119 129 293 49 86 102 270 31 164 64 52 72 71 78 285 317 50 108 63 79 45 70 97 148 126 69 53 100 107 163 304 117 83 97 314 329 146 102 115 131 62 65 56 146 281 113 134 60 245 161 352 211 45 150 67 51 62 102 117 138 99 130 440 40 18 29 100 156 134 165 98 54 41 212 77 98 161 121 89 72 124 107 208 125 235 92 261 280 48 53 209 254 102 197 18 113 128 208 335 162 103 105 64 244 76 305 209 149 93 203 98 195 220 69 147 93 240 136 309 30 85 36 44 46 82 162 194 372 64 204 67 120 110 137 97 163 203 339 69 67 199 47 98 278 218 194 203 106 57 194 324 385 74 469 111 64 58 78 120 345 67 242 61 191 113 265 97 682 122 34 101 110 52 122 137 42 63 93 79 41 107 68 351 254 61 222 125 170 399 293 22 58 545 142 180 327 201 152 111 756 212 531 47 52 135 158 171 189 50 84 148 90 64 175 268 74 406 61 161 223 197 149 169 87 88 130 45 47 203 457 55 222 45 32 46 254 126 207 425 106 182 140 101 366 102 125 28 161 414 178 113 148 180 84 128 168 82 147 212 59 50 76 73 126 87 100 27 94 261 165 131 71 391 239 159 278 182 149 121 82 105 117 260 179 105 174 183 159 67 87 225 113 84 33 77 37 197 159 181 34 74 241 83 119 160 253 57 135 257 323 213 105 122 121 292 60 124 96 240 160 88 124 229 100 96 114 75 28 275 67 97 68 58 81 614 187 196 88 91 146 173 60 101 215 93 188 35 305 215 158 36 328 230 138 120 163 207 141 95 31 91 148 350 181 95 18 65 48 422 376 102 155 115 98 181 150 61 104 114 212 154 403 312 100 70 135 34 71 73 48 321 172 78 92 43 100 325 464 133 136 128 247 193 93 244 198 51 30 114 288 144 136 59 110 260 24 90 51 59 137 24 129 213 101 143 354 64 138 263 50 222 46 163 94 107 53 170 89 198 12 45 60 54 201 198 77 43 109 76 340 215 141 568 435 175 68 115 61 124 28 227 125 111 55 107 152 149 35 84 57 313 241 76 126 79 49 83 146 214 329 66 323 107 52 101 29 165 101 123 102 189 127 156 132 72 71 319 198 240 93 406 110 117 149 107 52 38 410 112 397 93 32 144 105 103 142 119 205 426 62 288 84 21 67 132 30 65 316 184 445 119 103 107 73 132 165 99 141 60 225 74 74 316 121 267 73 377 68 62 34 64 122 230 164 197 55 204 180 120 97 87 161 95 269 98 374 103 178 124 91 103 182 265 60 62 280 92 93 34 303 108 62 122 117 150 228 306 129 86 207 30 86 70 68 97 68 22 72 233 162 178 202 103 67 103 210 174 161 112 132 671 82 84 123 104 465 109 96 187 41 71 141 195 158 145 96 107 127 160 201 291 200 258 96 109 90 162 56 60 55 131 31 135 59 110 161 232 56 163 193 107 154 299 269 273 26 100 194 83 41 63 39 143 182 573 92 115 234 494 79 57 167 166 166 47 310 235 202 50 88 133 155 52 131 262 360 312 70 203 507 425 349 127 277 56 132 411 65 79 297 294 120 304 88 113 161 119 114 99 102 94 217 61 180 154 238 157 229 225 65 135 154 93 130 101 142 96 44 110 129 64 353 54 72 57 144 148 56 84 95 81 207 128 178 48 73 174 78 309 825 49 144 89 149 33 218 190 173 746 89 337
48 5 11 3 9 4 7 7 16 28 1 15 16 13 6 1 11 2 3 10 6 2 1 7 16 9 9 3 7 7 21 12 12 9 3 7 10 10 16 7 10 6 12 3 12 10 12 6 4 12 17 14 15 2 13 4 6 5 13 5 6 8 15 17 27 8 5 2 9 16 1 7 1 22 24 46 7 11 7 6 4 4 25 12 3 4 10 1 3 19 12 19 2 6 3 6 8 12 1 1 23 10 2 10 7 9 4 6 4 9 18 21 15 9 19 45 2 9 7 13 7 20 8 22 2 6 11 3 1 2 4 14 2 0 11 12 16 4 16 15 2 8 4 10 10 5 34 13 5 7 2 0 8 5 4 7 13 6 2 2 2 2 7 17 19 9 28 11 1 8 2 1 7 13 2 9 8 9 6 11 10 15 9 3 9 6 9 9 5 14 11 2 18 35 4 14 19 6 6 5 4 17 9 4 9 13 30 5 7 5 8 8 5 6 5 1 39 4 6 11 3 16 4 13 10 9 5 3 14 11 6 26 6 6 4 9 4 4 11 2 3 6 13 5 10 8 13 8 8 3 2 2 9 10 14 5 6 5 7 3 28 1 3 48 2 3 45 9 4 10 2 10 8 14 7 8 7 4 4 7 5 14 30 2 8 2 9 13 7 7 5 8 8 13 4 5 10 11 12 26 2 6 2 13 2 11 7 1 15 2 7 2 0 28 9 2 26 5 4 4 2 8 5 6 8 12 7 13 12 22 8 14 11 4 7 28 6 6 3 4 6 7 15 11 2 5 8 9 5 6 14 26 6 13 4 1 8 4 4 2 19 10 7 14 24 8 4 5 5 17 6 9 6 19 4 3 5 2 3 11 5 17 15 4 7 7 16 20 13 3 12 3 4 4 14 8 5 12 2 13 9 6 15 3 5 23 18 13 21 4 8 7 3 17 8 16 15 8 5 16 12 17 11 5 5 6 4 9 7 19 2 3 3 1 12 21 4 14 3 14 29 8 3 3 1 14 5 29 2 8 6 8 11 17 19 7 12 30 5 0 12 21 3 5 14 6 6 5 13 7 19 2 9 2 3 14 8 8 5 6 7 6 9 13 9 9 3 2 8 45 10 3 6 1 4 3 6 5 8 3 20 4 7 9 15 3 17 5 22 7 4 11 11 3 2 12 2 7 3 7 3 15 4 9 7 7 20 38 19 3 7 2 3 13 15 10 23 16 4 8 9 5 3 8 6 4 9 6 14 9 25 7 3 27 7 5 9 4 14 3 8 2 1 0 1 1 3 2 9 8 3 7 14 9 3 8 5 9 11 8 6 7 2 9 9 4 2 10 14 11 9 22 10 13 0 4 8 7 9 5 15 20 10 9 14 11 13 1 20 5 5 9 6 2 6 16 4 3 8 5 1 24 5 10 21 7 2 12 5 1 10 5 2 5 5 5 5 4 13 0 5 5 13 11 6 22 14 1 9 6 10 9 20 3 6 4 10 2 12 24 5 5 19 30 15 12 32 5 2 15 1 12 5 5 2 13 2 6 2 9 25 7 8 9 2 37 13 22 5 0 10 5 36 8 11 2 11 7 9 4 13 1 5 5 5 8 3 19 19 5 5 5 10 1 6 15 18 4 5 8 10 4 12 22 3 23 7 6 16 4 4 9 7 6 3 16 2 17 10 5 6 7 22 16 9 9 12 11 25 8 10 10 8 11 7 9 3 12 5 18 2 10 18 9 7 4 14 9 9 17 6 12 16 7 7 4 1 7 0 14 5 1 8 6 7 17 15 13 19 1 0 19 8 11 22 4 3 4 3 8 12 14 9 3 33 2 5 5 4 18 2 4 20 7 4 31 11 6 12 21 11 13 6 15 3 34 12 4 3 13 2 18 4 13 9 6 8 7 4 9 10 4 6 4 8 4 14 5 3 10 14 12 5 13 16 7 8 29 11 2 7 14 5 9 2 18 3 19 12 1 7 11 2 5 7 5 10 5 4 7 13 12 13 2 3 10 18 14 13 5 6 2 12 2 25 5 1 10 4 26 7 4 9 9 1 9 3 3 2 12 17 24 10 6 4 9 33 15 24 20 6 8 5 4 8 3 19 4 5 15 37 7 14 21 15 9 6 10 4 13 6 3 14 11 21 5 14 4 4 9 4 4 11 7 2 7 5 1 19 3 14 12 32 13 12 3 28 10 15 7 4 14 31 2 2 5 8 13 18 9 31 28 1 12 1 8 5 11 17 6 12 12 7 4 8 8 28 4 18 17 10 4 10 5 26 12 34 11 13 25 3 8 6 1 3 2 14 3 8 37 14 18 3 5 4 4 17 9 8 9 25 13 6 4 0 14 9 14 2 7 11 3 6 4 2 5 6 7 12 12 12 14 15 15 16 9 3 1 5 9 4 5 12 1 0 8 15 7 17 4 7 11 2 4 24 14 6 3 5 1 4 10 2 2 35 5 1 10 3 14 8 28 7 7 6 2 5 6 15 6 6 9 7 11 10 14 7 5 9 6 14 12 8 22 8 14 20 36 13 9 5 1 5 23 3 9 8 4 10 5 7 8 28 25 5 12 9 4 2 5 12 18 9 5 23 14 6 5 4 29 4 5 14 30 15 2 4 2 14 14 13 2 12 11 1 5 0 3 10 12 7 17 6 10 3 11 2 26 13 8 10 2 8 23 10 6 7 1 10 8 1 9 0 4 2 11 18 7 21 10 4 1 14 37 11 21 7 5 4 6 31 11 5 49 4 12 19 3 4 12 10 6 12 4 13 13 1 7 7 1 18 13 4 9 13 20 0 13 5 15 7 2 7 19 15 16 23 11 7 25 13 12 6 9 11 20 5 7 4 24 11 12 7 4 12 4 8 6 10 9 16 3 3 7 10 13 7 14 6 4 5 17 16 4 6 1 9 6 4 4 6 1 9 1 11 11 17 13 1 4 8 16 5 3 3 6 26 7 12 7 6 5 7 5 6 2 10 9 7 16 10 5 6 7 14 6 4 15 3 6 7 4 29 11 1 5 16 1 6 9 4 5 6 5 2 6 14 5 15 16 3 8 14 20 1 7 9 14 6 2 18 16 11 10 1 17 12 14 1 11 17 19 1 6 6 6 8 14 2 10 9 5 9 17 2 25 15 5 17 14 23 16 20 4 3 7 13 11 5 8 13 9 22 21 8 4 5 7 22 1 5 4 4 33 12 6 4 9 10 9 2 14 1 12 6 4 3 4 5 10 4 6 3 3 3 9 9 9 5 26 3 2 8 2 4 18 1 20 5 4 9 11 16 3 3 7 7 2 10 5 3 14 2 30 4 3 4 10 12 10 15 5 8 16 0 2 17 10 3 10 2 7 2 6 5 3 11 17 3 22 12 4 8 10 12 5 10 1 5 8 9 14 5 13 4 9 13 5 9 25 0 24 11 17 13 20 9 3 1 2 1 5 6 24 5 5 9 9 20 12 25 22 5 5 11 4 2 9 4 12 4 40 14 8 12 19 9 6 8 7 2 2 13 6 7 3 1 8 8 6 13 1 3 5 3 10 4 7 11 2 9 7 7 9 27 21 2 2 6 1 7 1 8 12 4 26 16 1 6 7 5 3 12 2 1 6 15 2 26 4 5 12 10 10 10 29 8 3 18 5 4 1 26 7 4 8 8 9 4 15 32 12 19 8 2 5 5 4 12 15 15 20 3 3 10 12 5 7 2 4 4 5 6 10 24 9 9 4 6 5 19 4 4 12 4 7 7 16 4 23 9 3 7 3 13 3 1 8 6 2 4 5 7 14 2 6 8 38 7 3 4 4 3 7 13 7 9 3 14 61 7 17 7 24 15 3 10 17 8 1 12 15 8 20 3 20 10 7 13 6 8 7 7 2 9 1 10 3 9 10 3 3 7 6 14 4 13 7 15 18 6 13 39 8 4 3 14 8 27 3 3 9 19 14 5 8 4 6 3 13 46 6 7 2 21 5 30 7 9 14 4 10 8 3 9 10 3 9 4 25 3 8 15 3 30 20 24 2 27 6 21 4 3 7 6 10 33 1 3 4 11 12 29 11 1 1 12 3 3 7 6 11 3 16 3 3 7 10 16 5 23 10 15 10 4 11 3 4 10 8 11 7 2 5 7 4 8 6 15 2 9 9 8 21 30 16 60 19 10 8 16 2 28 4 7 1 27 10 13 6 19 17 10 22 7 71 4 12 11 6 4 15 9 6 11 2 14 4 20 18 16 1 9 5 3 8 6 13 5 9 18 26 16 8 11 8 18 15 12 17 11 3 7 8 36 10 11 18 4 5 19 10 19 16 10 16 8 6 8 6 6 4 5 5 7 10 6 4 13 16 5 10 14 12 9 5 6 19 13 4 2 5 11 0 16 6 6 7 5 2 8 8 2 10 4 25 7 22 2 5 2 3 8 6 2 8 3 2 13 1 17 10 8 9 0 4 16 6 6 7 7 6 6 9 7 9 5 4 3 7 2 7 2 13 9 9 6 0 6 18 4 3 1 4 17 10 5 9 16 8 3
90 12 29 20 8 15 13 9 23 17 9 12 13 11 5 13 19 9 27 3 6 36 22 8 13 22 9 11 8 46 28 36 2 29 7 11 2 7 25 15 4 7 31 32 44 12 25 8 11 7 18 13 14 2 18 8 8 13 6 8 7 19 53 23 5 13 13 6 13 49 4 15 20 12 20 7 8 6 15 21 2 2 14 10 27 3 19 15 11 6 10 28 40 12 15 21 6 24 4 17 29 24 8 45 17 20 26 8 11 7 15 15 6 12 10 22 8 12 32 34 16 17 11 18 43 17 24 11 4 15 4 12 11 6 19 19 24 8 16 20 9 3 12 20 16 7 14 20 3 5 8 10 27 32 13 12 14 17 7 15 5 6 8 12 30 16 16 18 32 8 9 6 22 34 9 23 3 18 23 10 19 15 25 19 6 4 3 20 3 7 8 9 5 36 1 20 19 22 25 26 16 20 29 11 82 32 5 27 51 18 4 9 32 2 26 20 11 30 16 14 12 14 12 4 6 25 45 10 17 11 19 17 8 11 17 18 20 9 11 10 13 6 17 9 2 9 12 13 7 20 14 14 9 82 12 36 13 9 9 10 21 16 7 12 56 53 11 7 7 8 9 7 20 7 30 20 9 8 11 21 19 8 43 7 15 13 8 2 9 20 4 23 7 24 20 6 21 19 17 10 12 28 1 22 9 14 4 6 10 14 21 8 6 34 39 19 6 12 4 11 22 8 25 52 9 5 10 21 23 7 54 3 13 11 8 16 24 22 5 15 10 27 10 7 24 15 13 49 7 7 7 10 28 5 12 8 5 31 11 22 18 20 7 16 33 21 7 16 22 4 4 9 17 5 12 21 18 19 20 9 32 11 12 20 6 16 4 11 32 14 16 17 26 22 18 13 5 31 9 15 12 8 7 17 10 12 41 9 4 14 18 15 11 16 20 24 9 3 4 35 11 8 3 86 7 4 3 10 5 16 17 17 8 8 52 9 8 18 19 7 15 8 28 8 26 21 17 6 18 6 33 3 47 9 12 12 15 36 6 12 9 6 12 15 6 54 14 27 12 9 10 12 6 5 12 21 9 45 12 14 6 17 20 25 9 19 4 21 10 39 4 9 17 5 12 9 9 10 14 10 32 16 12 25 49 11 11 13 21 7 16 42 23 9 7 14 20 12 20 15 26 17 18 31 39 24 17 18 3 23 11 8 4 5 19 17 25 10 16 3 14 11 16 12 3 6 6 7 35 6 11 14 10 21 13 14 17 8 5 17 5 12 3 7 10 6 33 10 3 9 9 3 26 10 8 4 8 14 31 10 18 4 34 21 22 8 18 8 32 16 10 13 9 15 15 12 4 6 13 7 9 23 13 8 35 15 0 20 5 16 7 30 6 24 21 14 5 8 9 34 18 25 23 8 14 9 26 3 21 18 12 10 5 21 25 88 8 3 6 30 15 11 11 2 5 46 12 7 12 10 29 14 9 22 9 18 11 1 18 14 11 15 13 35 14 91 11 19 29 11 69 13 9 7 32 11 34 3 3 16 28 27 45 26 14 68 21 17 45 18 19 8 20 5 10 17 18 6 4 10 16 3 7 10 3 15 49 9 21 12 16 4 3 18 24 6 16 6 19 4 33 15 11 19 24 21 6 9 7 7 16 32 13 14 13 50 6 25 10 12 6 5 25 31 14 22 3 13 13 12 9 15 7 5 18 13 15 9 6 11 12 9 21 6 3 14 6 19 9 14 8 26 4 9 9 39 18 7 8 8 34 19 4 9 13 11 15 24 33 9 14 26 4 66 10 3 31 5 8 13 5 13 23 7 26 10 20 17 9 22 15 71 21 10 8 11 38 13 16 27 37 24 20 13 44 26 9 18 9 8 6 27 14 13 76 15 17 9 14 12 5 9 32 26 18 17 14 12 8 12 13 25 30 4 12 17 4 45 19 47 8 19 7 11 9 13 6 5 45 5 27 7 5 6 16 7 2 8 17 20 19 16 18 4 11 12 17 20 11 15 22 4 18 12 28 11 28 16 7 13 15 27 8 11 40 45 10 8 10 12 69 24 19 9 16 23 8 12 11 41 21 20 14 20 23 23 42 11 10 8 34 12 10 6 21 26 7 8 43 10 8 14 13 19 13 7 6 10 38 5 11 13 10 8 17 8 11 5 29 5 7 6 8 3 6 3 27 27 25 12 0 24 1 21 9 5 28 6 6 37 28 5 54 14 21 8 12 10 26 10 27 12 8 11 5 0 6 10 9 22 5 3 9 18 11 2 10 4 6 22 10 8 29 16 17 7 11 9 3 15 7 14 23 17 27 6 4 2 23 37 5 33 16 24 40 30 6 22 10 2 2 32 12 4 16 9 10 46 8 11 17 9 4 18 8 54 28 36 10 51 6 6 15 15 4 15 12 29 27 10 13 13 3 13 13 13 10 26 4 9 13 6 19 20 16 40 5 14 25 20 49 11 6 13 30 5 38 18 16 11 20 4 42 7 17 18 5 1 5 4 18 17 12 8 5 21 10 10 11 7 25 39 23 15 3 17 12 16 13 19 14 20 4 8 10 18 15 4 9 8 11 15 16 51 23 7 5 13 26 16 2 17 27 14 11 8 9 44 16 22 18 45 7 58 25 16 16 14 10 8 5 12 35 9 22 6 24 13 7 32 13 32 0 10 4 13 11 10 1 19 40 14 18 13 3 5 8 28 6 24 20 40 12 6 28 23 13 11 10 23 18 9 2 15 8 9 10 15 12 17 15 11 14 12 26 16 7 8 20 1 18 39 9 14 9 4 2 19 7 24 11 4 53 6 5 8 39 12 6 12 31 13 6 16 6 40 5 12 11 11 5 4 12 4 17 10 8 1 21 12 11 21 25 22 58 43 9 23 26 7 13 11 6 12 18 22 54 10 8 17 28 29 12 15 3 31 7 13 15 50 10 9 20 8 4 19 3 21 20 3 18 30 26 17 33 31 35 5 24 15 18 3 6 14 11 3 25 8 22 2 30 9 5 12 11 57 24 5 49 13 15 83 28 21 12 8 18 37 14 27 25 8 25 8 12 17 11 13 15 8 2 7 30 9 16 16 15 13 13 53 5 3 17 6 28 10 18 18 14 6 6 7 6 8 34 37 25 20 43 35 13 55 42 17 20 12 33 11 11 34 11 1 39 7 45 2 30 6 13 16 11 20 29 24 28 15 21 8 12 32 5 33 19 33 7 31 10 19 15 15 9 7 20 14 31 32 10 3 9 25 14 13 22 16 5 8 11 33 6 31 13 24 27 0 8 5 2 8 9 23 15 27 11 18 43 33 14 9 16 5 7 19 6 10 4 38 9 19 8 25 8 9 26 30 16 15 24 14 22 34 13 15 11 21 11 2 16 2 9 7 25 6 24 14 4 21 25 55 11 4 16 9 15 4 12 7 9 5 9 17 23 9 17 13 46 2 17 9 6 9 11 16 15 5 14 10 22 25 42 21 4 30 17 21 15 8 33 3 3 26 9 2 8 9 7 11 4 7 33 14 11 8 14 28 29 6 8 5 13 4 16 30 13 22 24 17 5 10 9 34 22 6 23 17 4 8 17 12 20 7 9 6 7 21 10 7 22 15 4 6 10 16 33 26 21 20 14 9 14 11 18 19 16 24 7 49 18 67 13 74 18 8 33 6 11 17 21 15 15 10 25 12 8 17 29 7 18 3 11 14 19 12 40 1 39 9 7 17 6 21 10 52 10 6 12 9 100 13 14 7 6 12 41 15 31 10 5 13 32 9 17 2 11 5 66 20 18 9 4 5 22 27 9 9 12 19 22 15 16 18 11 14 10 9 10 5 19 8 15 14 53 10 42 30 14 9 26 19 7 17 8 15 11 7 13 8 3 18 19 13 39 6 11 0 4 2 13 17 4 17 25 14 15 9 13 11 1 6 18 9 19 38 11 3 10 16 11 10 12 16 12 13 1 23 13 18 11 8 23 3 22 11 6 5 6 8 5 21 27 15 2 13 5 6 39 5 27 10 20 6 2 7 15 12 33 5 7 2 11 14 29 9 27 15 2 5 11 7 16 13 5 21 5 23 33 20 5 39 15 5 5 9 66 5 9 25 68 8 4 24 19 19 16 4 9 7 10 6 22 6 8 25 11 16 22 7 23 7 22 7 18 5 4 4 8 11 12 36 11 19 24 11 25 11 12 27 11 20 10 13 4 9 12 5 2 14 18 40 4 21 8 15 9 14 14 25 6 10 17 16 12 9 25 13 10 13 44 8 8 3 15 9 9 10 17 7 38 9 3 12 22 8 32 24 12 7 11 25 27 6 15 17 10 18 9 20 8 13 6 44 17 7 24 21 9 12 30 11 63 45 13 9 4 19 23 35 4 3 3 12 11 20 13 8 39 39 65 28 16 8 14 22 16 10 12 63 26 4 37 18 6 7 19 15 5 11 13 3 13 27 2 26 4 16 30 8 9 7 10 8 11 18 8 20 33 24 34 4 28 8 30 9 25 14 17 15 11 18 2 6 3
91 4 5 1 6 15 2 4 3 6 8 9 14 6 2 5 8 3 4 3 1 0 6 10 2 9 1 7 6 13 7 10 25 8 11 5 7 1 5 2 2 2 1 2 5 7 2 12 17 7 4 3 5 5 2 1 1 10 15 4 5 4 15 5 17 3 8 3 7 8 2 5 9 0 23 5 7 4 11 9 2 8 5 0 4 2 13 0 6 3 5 5 1 5 7 3 17 19 7 12 4 4 3 12 10 6 11 11 0 0 4 4 28 10 6 2 13 7 8 2 4 2 2 9 11 13 47 12 13 1 5 7 3 0 7 4 6 9 2 2 7 4 9 6 6 1 2 17 0 4 9 5 3 2 16 6 2 1 14 6 1 3 0 7 10 1 32 15 14 13 10 28 5 19 12 7 4 3 13 9 2 0 5 31 9 4 16 5 8 9 2 3 2 16 4 15 2 4 17 2 8 6 6 9 0 11 7 2 10 7 6 3 3 7 2 3 8 4 49 10 4 13 2 10 5 6 3 2 3 4 12 4 11 8 6 3 2 6 10 23 6 7 15 2 8 4 4 11 10 11 7 10 4 0 4 15 19 1 4 4 3 6 13 7 4 18 4 2 1 2 1 2 4 8 3 8 1 3 4 5 7 3 3 6 2 2 2 10 10 4 10 2 0 6 5 9 9 5 9 12 0 6 4 0 8 7 5 6 11 14 4 6 3 16 8 15 19 5 7 3 8 5 2 10 3 4 3 2 10 3 3 7 3 6 4 8 2 10 10 15 9 12 1 7 1 8 8 15 2 3 2 10 14 4 2 7 0 13 8 7 4 7 5 12 4 8 6 28 0 3 5 5 5 16 6 6 5 3 2 4 3 7 4 6 5 13 7 9 8 5 10 34 5 7 3 3 1 1 2 8 4 1 3 1 2 7 29 9 5 1 5 6 10 2 4 11 3 13 10 7 43 13 3 1 2 3 15 8 6 13 9 2 3 14 8 7 3 2 8 6 15 16 3 4 1 10 1 3 3 3 0 5 3 5 0 28 10 10 2 1 2 5 14 5 6 13 2 1 19 2 4 24 8 9 3 7 15 4 1 3 5 1 6 15 11 0 5 4 7 5 4 2 5 1 4 15 6 1 7 9 23 9 22 12 3 9 3 3 10 8 7 10 4 8 12 3 2 4 4 8 8 7 8 15 6 50 14 8 2 6 25 2 0 4 1 9 3 3 4 4 4 7 5 4 3 6 3 4 9 18 2 1 6 11 10 4 6 19 0 2 3 7 5 14 3 1 26 1 1 8 7 16 9 6 1 4 4 5 1 2 6 15 9 8 6 1 6 2 11 6 1 4 10 2 4 6 17 2 20 2 5 2 16 10 15 8 2 5 4 3 4 3 4 6 10 2 12 11 5 14 2 7 15 2 7 2 7 16 11 5 3 3 16 11 5 2 4 7 4 6 7 10 18 8 5 11 3 10 2 5 2 12 2 2 10 10 1 4 10 6 7 4 10 4 6 4 1 5 11 2 5 13 20 37 1 4 7 2 3 5 14 4 10 23 2 9 2 8 6 15 8 2 3 5 5 11 5 1 9 4 1 1 4 9 9 8 5 6 3 2 3 4 19 19 15 17 14 6 7 11 4 4 4 6 4 2 20 2 9 7 7 3 10 9 1 5 16 21 5 3 3 0 4 12 12 8 6 2 4 20 3 10 6 9 2 9 5 2 25 8 4 11 5 1 2 16 5 2 7 7 13 2 4 10 19 8 20 3 3 5 9 2 7 4 8 5 3 1 5 14 7 6 15 6 12 17 2 14 5 6 3 7 4 17 5 3 19 8 5 6 7 5 10 1 11 8 3 4 9 11 10 4 6 3 11 17 8 5 4 11 4 3 5 6 8 20 8 1 15 10 3 7 7 13 4 15 5 8 1 1 3 7 14 14 20 4 13 4 6 10 1 3 8 9 10 6 16 11 11 15 4 5 8 3 2 13 24 5 1 11 6 8 3 1 4 1 10 10 3 7 7 3 8 4 9 3 2 8 3 14 2 5 10 10 26 6 3 9 0 4 0 10 2 22 12 3 2 6 11 6 13 2 10 10 7 9 8 5 14 8 2 7 6 10 6 5 0 3 2 10 7 12 5 3 8 7 4 5 4 0 1 2 14 5 2 8 22 5 3 4 10 5 10 3 2 3 35 3 3 6 5 2 1 7 2 15 8 3 13 2 5 10 12 2 12 0 5 8 9 7 7 4 5 4 6 5 9 2 5 23 7 6 12 11 4 6 5 13 7 2 3 5 3 23 11 13 0 1 0 6 1 7 4 4 10 0 5 5 6 7 5 6 17 5 3 4 5 1 0 1 0 4 11 4 7 8 3 1 0 2 8 12 2 4 4 7 3 15 18 12 7 2 10 14 4 7 1 2 6 11 2 12 11 4 1 6 9 2 2 6 6 4 5 12 4 2 3 4 11 4 7 11 1 5 12 12 13 5 4 5 2 3 0 4 10 13 12 5 13 15 2 21 12 9 9 2 0 7 6 3 11 3 6 13 3 9 8 10 11 4 18 1 4 8 49 10 8 13 6 9 14 4 9 0 20 5 4 15 8 16 0 17 5 8 7 3 4 4 12 4 6 16 11 6 6 3 5 3 3 8 5 6 3 3 0 13 7 11 1 5 3 4 1 1 10 2 5 11 12 5 13 3 7 10 13 0 6 18 4 24 3 5 11 12 1 7 13 10 14 1 7 8 2 4 14 11 11 3 7 6 6 5 0 3 11 2 3 2 5 3 6 9 5 6 3 8 0 4 2 5 8 1 33 4 6 9 5 6 11 5 4 9 11 1 5 1 2 9 9 11 2 2 3 9 10 7 0 9 1 3 3 8 8 5 12 4 0 17 7 4 6 11 4 3 18 3 1 14 10 7 6 2 8 6 2 11 3 11 13 9 6 2 11 2 9 0 10 2 3 7 7 22 5 1 5 14 4 8 2 10 7 1 4 14 8 3 7 13 2 9 10 4 12 14 5 9 6 3 9 2 7 0 2 7 4 4 17 1 6 6 14 6 11 3 13 4 5 2 16 4 5 3 41 4 6 10 8 7 4 3 3 19 1 8 1 6 4 0 6 1 6 1 5 7 1 15 6 7 7 4 9 5 7 8 1 8 20 5 5 12 3 17 3 5 7 6 3 10 2 5 8 3 0 10 12 8 12 8 5 7 4 16 5 20 11 1 2 9 2 13 11 3 7 7 34 1 10 4 1 7 7 9 11 3 4 4 9 12 3 4 4 13 13 1 27 9 4 6 11 15 0 3 8 9 1 5 10 11 1 11 4 7 9 7 1 15 12 2 1 3 4 1 6 8 4 2 7 11 17 0 5 3 2 15 3 2 3 2 0 2 3 23 13 2 3 2 4 4 6 5 3 2 1 14 3 3 6 3 9 3 6 4 3 16 33 4 1 2 7 2 5 10 1 18 6 3 9 2 7 6 10 0 14 14 2 20 6 4 5 1 17 18 6 14 1 5 0 0 8 1 12 2 11 23 14 7 15 2 7 6 1 15 4 4 21 6 5 10 6 13 2 6 4 5 5 11 8 0 6 11 9 1 3 2 3 1 2 7 4 11 6 3 21 10 6 8 5 6 3 8 9 45 9 4 8 3 2 5 11 10 1 3 1 2 4 9 5 2 4 10 13 8 4 4 2 5 38 5 2 5 8 5 6 4 27 10 8 2 6 13 8 10 0 3 3 0 8 10 17 7 2 1 2 8 8 4 8 7 10 1 5 13 9 6 6 1 3 3 5 3 1 11 2 4 1 2 1 4 6 3 12 6 5 7 23 9 15 14 10 18 5 16 17 15 17 4 5 8 3 20 12 14 0 3 2 3 2 4 6 5 5 5 2 1 7 0 3 4 12 7 8 15 5 8 6 5 21 3 0 7 4 13 3 5 5 2 12 0 4 6 5 0 14 5 8 12 8 12 11 0 1 13 8 9 10 8 8 1 6 7 14 6 9 3 2 12 12 8 5 3 12 3 7 9 4 11 5 1 34 8 3 5 13 2 2 1 6 5 2 7 0 16 12 9 15 10 2 1 10 16 2 14 21 3 6 4 6 0 1 5 6 12 3 8 5 8 3 1 4 10 2 26 3 13 8 10 6 8 1 8 26 10 16 3 6 2 10 14 21 7 7 17 7 5 15 9 9 4 6 1 6 7 10 3 9 0 11 6 4 4 3 4 4 7 9 20 8 12 5 4 10 3 1 3 4 11 11 20 8 7 8 5 5 5 10 6 4 9 4 1 6 4 6 3 10 3 8 16 24 1 3 1 5 4 0 1 17 3 19 20 17 2 11 16 7 2 2 4 2 5 8 10 5 4 18 5 8 8 4 5 12 4 1 3 4 3 6 11 7 4 11 4 5 14 17 11 10 2 4 15 3 0 6 4 6 3 0 2 9 17 4 6 13 1 2 10 2 2 1 7 8 18 16 12 4 21 8 5 9 3 4 4 10 2 4 7 10 4 2 7 1
208 7 4 8 1 4 8 2 12 2 6 5 1 24 3 10 12 1 3 2 1 12 3 4 13 4 3 4 7 2 6 7 5 6 10 4 4 9 6 9 2 7 1 9 15 4 4 7 11 2 8 4 4 2 8 3 0 2 6 1 3 3 5 9 9 3 1 4 5 0 3 8 2 6 5 5 4 3 7 7 4 6 22 4 6 4 0 0 0 3 1 7 0 18 6 5 7 4 5 0 4 6 3 1 4 18 3 3 9 12 13 2 0 0 3 3 1 2 3 7 1 4 8 8 11 7 2 7 4 6 0 3 8 9 1 3 5 4 10 1 1 12 3 4 1 3 6 3 4 2 4 1 11 2 13 1 3 6 1 4 5 1 5 10 3 2 3 4 3 5 9 3 4 6 7 7 3 3 13 6 2 1 10 7 2 1 2 7 2 5 4 4 5 23 2 1 2 4 7 8 2 4 2 1 1 1 2 2 4 9 6 5 5 2 1 10 2 6 0 3 2 8 0 1 4 4 2 6 5 13 0 4 1 8 7 1 4 4 2 9 5 6 5 1 13 3 2 3 7 6 4 4 7 2 4 0 9 4 9 3 2 3 3 16 13 3 0 3 5 13 3 1 4 18 4 3 30 5 2 5 2 4 6 3 4 6 8 4 2 1 6 5 6 4 3 0 1 5 2 1 0 7 1 15 1 1 4 2 6 2 9 3 6 16 6 3 4 4 3 1 1 1 3 3 3 2 1 1 10 4 3 3 3 13 0 5 0 3 2 7 2 8 1 2 4 4 1 3 8 2 3 3 4 4 2 2 4 3 5 2 4 2 0 6 22 1 3 2 3 6 1 2 3 8 3 6 3 19 1 2 1 4 4 4 1 2 1 10 4 3 11 9 3 14 12 1 2 2 1 5 3 2 7 4 0 1 0 5 4 0 4 3 3 2 4 8 3 5 8 2 4 4 16 16 2 36 1 3 2 12 2 8 7 2 1 6 1 4 1 2 0 8 3 4 5 5 0 6 4 9 3 4 1 2 2 3 16 2 5 1 6 16 4 1 6 5 4 13 7 5 6 13 2 5 0 4 1 8 12 4 21 6 3 7 3 3 1 2 7 8 4 5 7 0 4 6 8 2 12 0 43 17 5 5 8 3 2 5 1 14 4 5 4 10 2 1 1 2 2 6 10 11 0 23 16 2 12 6 4 2 5 1 1 13 0 11 2 4 5 2 2 12 8 3 8 12 5 1 11 7 6 9 3 10 6 3 0 1 4 2 1 23 2 1 0 6 3 8 7 1 3 1 4 8 3 1 5 1 3 3 3 1 6 2 4 12 4 8 7 4 0 3 4 3 7 3 5 6 3 1 8 1 3 3 13 5 17 3 8 3 0 9 7 4 11 5 15 12 2 2 2 8 3 2 7 12 2 13 1 1 6 9 1 2 4 2 7 4 0 4 1 3 6 32 19 18 2 1 3 1 14 14 4 3 6 1 3 1 3 4 0 1 4 3 8 13 5 4 3 7 5 4 2 21 7 6 1 5 2 1 3 6 7 8 1 4 4 5 4 1 4 3 5 3 0 4 5 1 2 2 3 1 12 2 1 2 3 17 17 0 4 6 13 1 2 3 3 2 0 7 7 8 6 3 5 2 7 1 8 12 13 7 0 6 10 1 5 6 1 24 1 4 3 6 1 3 10 5 3 2 1 1 10 2 4 3 4 4 5 2 1 4 7 1 12 3 7 11 3 5 1 12 8 3 17 1 3 0 6 5 8 1 5 1 0 7 0 2 5 1 1 4 3 4 4 1 2 25 5 1 3 3 14 0 6 14 3 3 6 3 3 6 10 2 0 3 4 2 10 4 9 11 0 3 2 7 4 9 16 7 6 0 6 7 1 5 2 1 1 0 7 6 1 3 5 1 3 9 3 5 2 4 10 3 67 2 11 3 4 5 1 2 2 8 2 7 3 2 3 2 9 8 3 3 0 2 8 0 5 19 9 8 1 4 3 14 13 5 1 2 1 7 4 1 9 2 4 0 1 10 5 2 5 8 21 7 2 0 7 2 4 3 2 9 5 3 2 7 1 7 0 8 8 7 5 2 3 1 9 17 4 2 5 6 1 4 4 6 4 4 4 5 1 7 4 0 6 4 11 5 1 7 5 3 4 5 3 1 6 0 1 6 4 2 1 1 0 1 4 1 1 3 3 4 6 5 1 9 1 2 1 1 0 5 5 3 3 6 7 5 2 3 3 3 3 2 0 7 2 2 0 3 6 11 4 3 5 6 2 2 4 2 6 8 5 5 4 4 0 2 9 14 7 5 6 1 5 0 3 5 4 7 1 13 3 10 7 2 3 4 6 1 3 2 7 5 2 7 4 9 3 0 2 4 3 0 2 8 1 1 3 4 1 4 9 4 16 11 4 3 5 3 4 5 4 4 3 6 2 1 2 5 8 1 6 2 1 1 4 2 6 10 7 2 9 3 3 17 3 3 2 3 3 6 4 1 7 1 22 0 2 5 7 2 5 8 12 4 3 5 5 4 11 4 6 6 2 5 6 0 2 5 14 8 5 2 10 2 2 2 11 3 3 4 3 7 2 6 3 7 1 5 9 8 4 10 12 2 13 8 11 2 0 9 1 1 0 1 2 3 1 2 5 2 3 5 4 0 3 3 5 4 9 5 1 0 4 3 2 1 1 8 6 7 8 6 7 4 5 5 4 2 4 1 7 8 3 5 4 8 39 0 6 8 4 9 1 7 3 4 2 6 2 1 4 13 6 5 3 0 6 7 8 23 4 4 3 1 2 5 5 5 6 2 2 6 9 6 1 4 1 10 0 1 6 3 2 7 2 4 3 1 0 1 1 12 5 2 2 6 14 3 4 2 10 19 3 2 3 3 5 3 7 2 2 5 6 7 14 3 2 1 11 5 3 3 7 16 2 1 1 18 4 6 2 2 1 3 7 3 4 3 13 1 12 4 2 4 3 4 8 11 5 1 3 1 16 3 6 0 2 0 3 8 2 8 8 4 7 2 4 4 5 8 2 4 18 7 1 2 1 4 8 1 6 6 6 5 1 8 13 3 4 5 1 5 3 5 4 3 4 11 9 2 3 7 2 11 6 5 4 0 4 3 3 1 6 3 13 3 3 6 5 12 2 0 4 8 6 5 2 5 2 4 2 3 8 22 2 9 2 4 1 4 10 1 1 12 7 5 2 6 5 1 1 9 2 19 1 8 3 4 9 6 6 7 5 9 13 6 8 2 3 0 6 5 21 4 1 2 8 7 17 6 1 6 1 2 20 5 0 0 2 3 1 6 2 2 16 8 2 2 3 1 3 1 12 7 9 0 1 0 1 5 0 2 2 2 11 1 2 12 2 3 4 1 3 5 1 0 1 1 2 0 3 2 1 2 4 4 8 2 5 2 5 3 8 6 6 9 4 0 1 6 7 19 4 11 4 3 4 7 1 2 7 4 2 7 1 8 3 3 3 8 13 8 1 1 5 4 9 5 2 7 1 0 5 3 7 1 7 7 7 6 3 9 2 3 6 1 9 5 2 1 2 11 4 0 4 3 0 8 2 7 5 5 2 3 5 6 6 1 3 14 1 2 1 2 2 2 5 8 0 6 11 1 4 9 5 3 4 1 5 6 3 4 12 1 2 3 11 1 4 3 6 3 1 6 6 4 3 2 1 2 9 10 2 9 6 1 4 2 3 4 2 13 1 2 1 19 14 5 5 5 6 4 3 9 1 4 2 5 2 1 2 8 1 5 5 0 1 5 5 5 5 2 9 15 3 5 0 0 3 2 8 2 1 9 7 5 4 2 2 4 4 6 4 3 3 5 1 3 15 19 2 10 2 9 3 2 6 2 3 4 3 0 0 3 13 2 3 4 4 4 4 3 7 12 4 5 1 3 2 2 4 3 9 8 10 4 2 5 7 5 5 0 2 11 2 6 4 1 3 2 8 4 4 7 2 14 6 1 3 2 7 2 9 4 7 6 10 7 8 6 2 14 9 2 7 5 5 6 1 1 3 3 8 8 2 6 5 5 7 1 1 8 3 2 8 8 2 4 3 5 12 13 6 8 6 6 0 10 9 6 2 2 6 1 5 2 6 11 10 6 4 5 6 1 8 2 1 4 10 3 1 4 2 6 9 12 5 2 9 12 1 4 4 7 2 4 14 1 8 0 33 1 2 4 3 20 1 2 3 4 20 3 2 2 5 2 2 5 11 2 8 2 2 10 3 9 5 3 2 11 5 20 1 7 0 3 8 1 3 5 4 4 0 0 1 4 3 2 7 11 4 4 13 5 4 2 3 5 8 0 5 3 3 4 6 5 3 5 0 8 3 3 7 6 3 14 1 2 4 8 4 8 4 4 7 6 0 9 1 24 2 7 3 3 6 2 6 2 2 2 5 13 5 8 4 9 0 4 4 7 3 11 4 6 10 0 10 5 2 4 4 5 5 5 5 1 6 6 0 3 5 2 2 2 4 10 9 2 2 0 6 9 5 2 4 3 1

Nota: La predicción de modelos requiere mucha memoria. Por lo tanto, solemos ejecutarla en paralelo para bloques de celdas de la cuadrícula.

Predicción de distribuciones

La modelización bayesiana es una modelización estocástica. De este modo, tiene en cuenta las incógnitas, causadas, por ejemplo, por:

  • Indicadores débiles

  • Errores en los datos de entrada

  • Representatividad incompleta de la muestra

  • Variación de las observaciones

Estas incógnitas generan incertidumbre en las predicciones. Más concretamente, permite predecir una distribución del número probable de habitantes.

Como se ha visto en la sección anterior, tenemos para cada celda de la cuadrícula 2000 predicciones. La imagen 2 muestra un ejemplo para cinco celdas de cuadrícula.

prediction_ex <- predictions %>% slice(1:5) %>% 
         pivot_longer(-gridcell_id, names_to = 'iterations', values_to = 'prediction') %>% 
  group_by(gridcell_id) %>% 
  mutate(
    mean_pop = paste0(gridcell_id, ' (mean=', round(mean(prediction)), ')')
  )


ggplotly(ggplot(prediction_ex,   aes(x=prediction, color=mean_pop)) +
  geom_density()+
  theme_minimal()+
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())+
  labs(x='Population predicted', color='\n\nGridcell id \n(mean population)')) %>%
  layout(margin=list(t=100))

Imagen 2: Ejemplo de recuento de población previsto para 5 cuadrículas y su recuento de población medio

Una estadística significativa para recuperar de la distribución es la media, que puede considerarse como el valor más probable.

Cuadrícula de población

La población ascendente en cuadrícula que publicó WorldPop se basa en el valor medio por celda de cuadrícula.

Veamos cómo podemos crear una población en cuadrícula a partir de nuestras predicciones.

Primero calculamos la predicción media para cada celda de la cuadrícula:

mean_prediction <- tibble( 
  mean_prediction = apply(predictions %>% dplyr::select(-gridcell_id), 1, mean)
  )

A continuación, añadimos la celda de cuadrícula no establecida uniéndola con raster_df, que contiene todas las celdas de cuadrícula.

mean_prediction <- mean_prediction %>% 
  # add gridcell_id
  mutate(gridcell_id = predictions$gridcell_id) %>% 
  # join unsettled grid cell
  right_join(raster_df %>% 
               dplyr::select(gridcell_id)) %>% 
  # sort according to position in the raster
  arrange(gridcell_id)

Para recuperar la estructura ráster, cargamos un patrón de referencia, normalmente el mastergrid.

wd_path <- 'Z:/Projects/WP517763_GRID3/Working/NGA/ed/bottom-up-tutorial/'
mastergrid <- raster(paste0(wd_path, 'mastergrid.tif'))

Y asignamos el recuento de población previsto.

# create raster
mean_r <- mastergrid
# assign mean prediction 
mean_r[] <- mean_prediction$mean_prediction

Por último, podemos trazar la trama para ver nuestra cuadrícula de población:

library(RColorBrewer)
# create color ramp
cuts <- quantile(mean_prediction$mean_prediction,
                 p=seq(from=0, to=1, length.out=7), na.rm=T)
col <-  brewer.pal(n = 7, name = "YlOrRd")

#plot mean prediction
plot(mean_r, col=col)
Gridded mean population prediction for region 8

Imagen 3: Predicción de población media en cuadrícula para la región 8

La imagen 3 muestra la distribución espacial típica de la población. Vemos en la esquina superior izquierda un grupo de ciudades, con barrios más densos en los centros. Las líneas de asentamientos más densos están relacionadas con la distribución espacial de la población en torno a las carreteras. En blanco se muestran las celdas de la cuadrícula no pobladas.

Incertidumbre en cuadrícula

A partir de la distribución a posteriori completa, también podemos obtener el intervalo de credibilidad (IC) del 95 % para cada predicción y la incertidumbre relativa correspondiente calculada de la siguiente manera:

\[\begin{equation} uncertainty = \frac{upper\_CI - lower\_CI}{mean\_prediction} \end{equation}\]

La incertidumbre relativa nos informa sobre las áreas en las que la predicción media es menos fiable.

Para calcular la incertidumbre en cuadrícula, seguimos el mismo proceso que para la predicción de la media en cuadrícula.

# retrieve the 95% credible interval
ci_prediction <- apply(predictions %>% dplyr::select(-gridcell_id),1, quantile, p=c(0.025, 0.975))

ci_prediction <- as_tibble(t(ci_prediction)) %>% 
    # add gridcell_id
  mutate(gridcell_id = predictions$gridcell_id) %>% 
  # join unsettled grid cell
  right_join(raster_df %>% 
               dplyr::select(gridcell_id)) %>% 
  # sort according to position in the raster
  arrange(gridcell_id) %>% 
  mutate(
    mean_prediction = mean_prediction$mean_prediction,
    uncertainty = (`97.5%` - `2.5%`)/ mean_prediction
  )

# create uncertainty raster
uncertainty_r <- mastergrid
uncertainty_r[] <- ci_prediction$uncertainty

#plot uncertainty raster
cuts <- quantile(ci_prediction$uncertainty,
                 p=seq(from=0, to=1, length.out=7), na.rm=T)
col <-  brewer.pal(n = 7, name = "Blues")
plot(uncertainty_r, col=col)

Se observa una mayor incertidumbre en las afueras de los asentamientos más densos. Esto suele ocurrir debido a la rápida expansión de la extensión de los asentamientos, que es difícil de registrar.

Predicción agregada

Las cuadrículas de población son ideales para visualizar la distribución espacial de la población en una zona. Pero, a menudo, queremos obtener agregados, es decir, una estimación de la población de una zona personalizada, como una ciudad o la zona de influencia de un hospital.

Veamos un ejemplo. Dibujamos manualmente la extensión de una ciudad de nuestra zona de estudio.

Advertencia: las extensiones no son las oficiales

library(sf)
library(tmap)
tmap_options(check.and.fix = TRUE)
city <- st_read(paste0(wd_path, 'study_city.gpkg'), quiet=T)

names(mean_r) <- 'Gridded population'
tmap_mode('view')
tm_shape(city)+
  tm_borders()+
  tm_shape(mean_r)+
  tm_raster()+
  tm_basemap('Esri.WorldGrayCanvas')

Obtener la predicción de la población media de la ciudad corresponde a calcular la estadística zonal de la cuadrícula de población, es decir, sumar todas las celdas de la cuadrícula contenidas en la ciudad.

mean_city <- extract(mean_r, city, fun=sum, na.rm=TRUE, df=TRUE)
mean_city %>%  
  mutate(features = 'city', .after=ID) %>% 
  rename("Mean population prediction"=Gridded.population) %>% 
  kbl(caption='Mean population prediction for the city computed with the gridded population') %>% kable_minimal(full_width = F)
Table 2: Mean population prediction for the city computed with the gridded population
ID features Mean population prediction
1 city 376610.6

Obtener la incertidumbre de esa estimación es más complicado. De hecho, no podemos sumar la incertidumbre de la celda de la cuadrícula para obtener la incertidumbre del agregado porque los intervalos de credibilidad se basan en el cálculo de cuantiles y los cuantiles no pueden sumarse.

Por lo tanto, necesitamos reconstruir la distribución de la predicción de población completa para la ciudad agregando todas las celdas de la distribución de la predicción de población.

Primero convertimos el polígono de la ciudad en un ráster con la misma extensión que la cuadrícula de población. Este ráster se compone de 1s para las celdas de la cuadrícula dentro de la extensión de la ciudad y de 0s para las celdas de la cuadrícula fuera de la ciudad.

city_r <- rasterize(city, mean_r)

Convertimos el ráster de la ciudad en un conjunto y lo unimos a raster_df para obtener el identificador de la celda de la cuadrícula. Filtramos las celdas de la cuadrícula que pertenecen a la ciudad y las combinamos con las predicciones:

city_prediction <- raster_df %>% 
  # select grid cell id from raster df
  dplyr::select(gridcell_id) %>% 
  # add city dummy
  mutate(city = city_r[]) %>% 
  # keep grid cells inside the city
  filter(city==1) %>% 
  # join predictions
  left_join(predictions) %>% 
  # keep predictions
  dplyr::select(starts_with('V'))

city_prediction contiene todas las celdas de la cuadrícula comprendidas en la ciudad con la correspondiente distribución de la predicción de población completa.

Primero sumamos las predicciones de cada celda de la cuadrícula en cada iteración

city_prediction <- as_tibble(apply(city_prediction,2, sum, na.rm=T))

city_prediction se convierte en una matriz de tamaño 2000 que contiene, por lo tanto, los totales de población de 2000 de la ciudad.

A partir de la distribución de la población total de la ciudad, podemos obtener estadísticas significativas.

city_prediction_stats <- city_prediction %>% 
  summarise( 
    "Mean population" = round(mean(value)),
    "Upper bound" = round(quantile(value, p=0.975)),
    'Lower bound' = round(quantile(value, p=0.025)))

ggplot(city_prediction, aes(x=value))+
  geom_density(size=1, color='orange')+
  theme_minimal()+
    theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())+
  labs(y='', x='Predicted population for the city')+
  annotate("segment",x=city_prediction_stats$`Mean population`,
           xend=city_prediction_stats$`Mean population`, y=0, yend=0.000005, size=1)+
  annotate('text',x=city_prediction_stats$`Mean population`+5000, y=0.0000012, hjust = 0, fontface =2,
           label=paste0('Mean prediction: \n', city_prediction_stats$`Mean population`, ' people'))+
  annotate("segment",x=city_prediction_stats$`Upper bound`,
           xend=city_prediction_stats$`Upper bound`, y=0, yend=0.000005)+
  annotate('text',x=city_prediction_stats$`Upper bound`+5000, y=0.000001, hjust = 0,
           label=paste0('97.5% prediction \nbound: \n', city_prediction_stats$`Upper bound`, ' people'))+
    annotate("segment",x=city_prediction_stats$`Lower bound`,
           xend=city_prediction_stats$`Lower bound`, y=0, yend=0.000005)+
    annotate('text',x=city_prediction_stats$`Lower bound`+5000, y=0.000001, hjust = 0,
           label=paste0('2.5% prediction \nbound: \n', city_prediction_stats$`Lower bound`, ' people'))
Full distribution of the predicted population total for the city

Imagen 4: Distribución completa de la población total prevista para la ciudad

La imagen 4 muestra que la predicción media coincide con la calculada con la cuadrícula de población. Además, vemos que nuestro modelo produce predicciones con un intervalo de credibilidad muy amplio.

Disponer de la distribución completa permite responder a preguntas como “¿cuál es la probabilidad de que haya más de xx personas en la zona?” o “con un 95 % de certeza, ¿cuál es el número máximo de personas en esa zona?”

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

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.
Laboratorio Nacional Oak Ridge. 2018. “LandScan HD: Nigeria.”