5 Bottom-up Models
Bottom-up population modelling methods (Wardrop et al. 2018a, Leasure et al. 2020i, Boo et al. 2022a, Darin et al. 2022a) use geolocated household survey data from a sample of locations to fit statistical models that estimate population sizes for unsampled areas based on relationships with spatial covariates. WorldPop develops customized statistical models for individual countries to make the best use of available survey data and to provide robust estimates of uncertainty.
This is a good approach when there has not been a recent or complete national census but there are recent geolocated household survey data available. This approach provides Bayesian estimates of uncertainty but requires more detailed input data and more time to develop the models.
5.1 Input Data
Bottom-up methods require a few key types of input data:
- Population data
- Settlement map
- Geospatial covariates
- Administrative boundaries
5.1.1 Population Data
Population data used in bottom-up methods generally should include counts of people in clearly defined georeferenced areas. A polygon shapefile with the boundary of each enumeration area and the total population within each area is ideal. There are a few potential sources for these data:
- Partial census results
- Microcensus surveys designed for population modelling (a random sample of locations where enumeration is carried out)
- Pre-survey listing data from routine household surveys (e.g. DHS, LSMS, MICS)
Point locations of buildings and/or households within enumeration areas are sometimes collected during census and survey field work. These data can be very useful because they provide higher resolution information about population patterns, although they are not required. Pre-survey listing data can also be very useful, especially if the surveys were recently conducted in areas that were inaccessible to census enumerators. If pre-survey listing data from household surveys are used, additional information about the site selection is also required. For example, If the household survey used a sampling design in which survey locations were selected with probabilities proportional to population size (PPS), then it will be necessary to obtain the weights used for PPS sample design.
5.1.2 Settlement Map
A settlement map identifies areas where residential structures exist. It may also classify areas into settlement types, such as urban, peri-urban, rural, slums, commercial and industrial (see chapter 4). This information may be in the form of:
- Building locations (points)
- Building footprints (polygons)
- Gridded map identifying pixels that contain buildings (raster)
These data could be derived from several sources:
- Satellite imagery
- Pre-census cartography
- Building points and footprints
- Gridded derivatives of building footprints (freely available for some countries (Dooley et al. 2020a))
If there is no classification of settlement types available, building points or building footprints can be directly used to identify different settlement types based on the patterns of building locations (Jochem et al. 2018, 2020). There are also freely available global settlement maps, but the quality from global data sets varies strongly between countries. The smallest settlements are often missed, hence this would need to be considered before committing to any publicly available global settlement map.
Additional data for each building such as building area, height, or use (i.e. residential, commercial, mixed) can be very beneficial for population modelling. Classifying individual buildings as residential or non-residential (Sturrock et al. 2018, Lloyd et al. 2020) can sometimes be accomplished with existing public data from Open Street Maps and other sources. While these additional data can improve population estimates, they are not required.
5.1.3 Geospatial Covariates
Geospatial covariates are spatial datasets with national coverage that describe any variable that may be correlated with population densities. There are many suitable datasets that are publicly available, including some produced by WorldPop for this purpose (Lloyd et al. 2019, Dooley et al. 2020a).
For example, a digital map of road networks (a line shapefile) could be used to calculate road densities which may correlate with population densities. Or, global satellite-derived nighttime lights data sets (raster files) may correlate with population densities in some areas. Administrative records could also be useful such as electricity usage for each administrative unit (polygon shapefile). Locations of public facilities such as schools (a point shapefile) can also be very informative. If the number of students attending each school is known, that would also likely add to the accuracy of population estimates.
There are an almost infinite number of possible geospatial covariates. Many of them are publicly available for population modelling. However, it is essential to identify good quality covariates (i.e. those that are strongly correlated to population density) and that has a comprehensive national coverage to significantly improve the accuracy of population estimates.
5.1.4 Administrative Boundaries
Administrative boundaries may include regions, states (provinces), divisions and sub-divisions. These administrative units are often nested within one another. Administrative units can be used by the model as a covariate to improve estimates of population densities. Administrative units can also be used to summarize model results, providing population totals for each administrative unit.
5.2 Statistical Models
WorldPop develops customized Bayesian models to make the best use of available data for specific countries and to accurately quantify uncertainty associated with the population estimates.
Bayesian models generate population estimates as probability distributions known as “posteriors”. Examples of posterior probability distributions for population estimates can be found within the woprVision web application. The mean value of the posterior probability distribution is taken as the expected value for the population estimate. Variance around the mean represents uncertainty in the population estimate.
Uncertainty in population estimates may be caused by several factors. Sometimes, uncertainty results from sampling error associated with small sample sizes (i.e. not many household survey clusters in the area). Uncertainty may also represent true variation in population densities from neighborhood to neighborhood that simply cannot be explained by the covariates in the model. Uncertainty may also relate to the structure of the statistical model itself. Uncertainty in a model can be reduced by A) collecting more household survey data, B) finding better covariates to predict population densities, or C) revising the model structure.This require a cost-benefit anaysis. Revising the model structure is by-far the easiest and is one reason why the flexibility of Bayesian models is so important.
5.2.1 Software
Before exploring the models themselves, it is important to be familiar with the required software. The R programming language for statistical computing (R Core Team 2020a) is ideal for fitting Bayesian models. There are a number of software packages available, but a few that WorldPop regularly uses are:
- STAN software (Carpenter et al. 2017) with the rstan R package (Stan Development Team 2020)
- JAGS software (Plummer 2003a) with the runjags R package (Denwood 2016a)
- INLA R package (Lindgren & Rue 2015)
For those new to Bayesian modelling, it is recommended to start with STAN. This is because it provides full flexibility to customize the models, it has excellent documentation (mc-stan.org), and it is computationally more efficient than JAGS. If users are already familiar with the BUGS or JAGS languages, all of the models described below can be built using either software. To build geostatistical models (see Geostatistical Models), INLA (R-INLA.org) is the recommended software, given that it is more computationally efficient than JAGS or STAN for estimating high-dimensional spatial covariance parameters. However, it is less flexible for building customized hierarchical models.
5.2.2 Simple Model to Start
A simple linear regression can be written as:
\[ y_i \sim Normal(\mu_i, \sigma) \\ \mu_i = \alpha + \beta x_i \]
where \(y_i\) is the value of the response variable at location \(i\) and \(x_i\) is the predictor variable (e.g covariate). These two variables represent the observed data and the rest of variables represent model parameters that will be estimated using STAN, JAGS, or INLA (software described above).
\(\mu_i\) is the expected value of the response variable based on the covariate value \(x_i\) at a given location. It is the mean of the normal distribution. Random noise that could result in the observed value being different than the expected value (for example, residual variance or uncertainty) is represented by \(\sigma\) (i.e. standard deviation). The regression coeffecient \(\beta\) (i.e. regression slope) estimates the effect of the covariate on the expected value, while \(\alpha\) (i.e. regression intercept) is the expected value of the response variable when the covariate is equal to zero.
The first line of the model is the stochastic model (i.e. it includes random noise) and the second line is deterministic (i.e. it always generates the same output for a given input). The selection of a normal distribution (a.k.a. Gaussian) in the stochastic portion of the model should be based on characteristics of the response variable. A normal distribution represents continuous numbers that can be negative or positive.
In population modelling, the response variables are the counts of people \(N_i\) which are always positive integers, as such, a more appropriate stochastic model is required. The above Gaussian linear regression can be modified into a Poisson regression as follows:
\[\begin{equation} N_i \sim Poisson( \mu_i ) \\ log(\mu_i) = \alpha + \beta x_i \tag{5.1} \end{equation}\]
This is a generalized linear model with a log-link function (McCullagh & Nelder 1989). The log-link function ensures that \(\mu_i\) is always positive, and the Poisson distribution produces positive integers. Now we have an appropriate deterministic regression and stochastic model for population counts.
5.2.3 Bayesian Priors
To implement this model Eq. (5.1) in a Bayesian context, priors must be defined for \(\alpha\) and \(\beta\). Priors are probability distributions that represent our prior knowledge about the range of possible values for parameters included within a the model. Priors must be specified for any “root node” parameters, those that do not show up on the left side of any probability statements in the model. Probability distributions used as priors are usually very disperse flat priors so that they do not influence the posterior parameter estimates. In general, priors that are informative enough should be specified. This ensures that a realistic range of possible values for the parameter are defined whilst remaining vague enough to allow for the observed data to have a dominating influence on the parameter estimates.
For the model in Eq. (5.1), uninformative flat priors can be used: \[ \alpha \sim Uniform(-10, 10) \\ \beta \sim Uniform(-10, 10) \] On the log-scale, this is a range from near zero to over 22,000.
Alternatively, more informative priors can be used:
\[ \alpha \sim Normal(0, 5) \\ \beta \sim Normal(0, 1) \] The relative influence of priors depends on the scale of the response variable and the structure of the model. It is good practice to test the relative influence of various priors on the posterior parameter estimates before making a decision on the prior to use.
Going forward in this chapter, priors will not be explicitly specified in the given examples unless the prior selection is noteworthy.
5.2.4 Hierarchical Core Model
A hierarchical model is one where the output (left side of equation) from one stochastic model serves as the input (right side) to another. Building on the Poisson regression in Eq. (5.1), we can build a hierarchical model that incorporates population density \(D_i\):
\[\begin{equation} N_i \sim Poisson( D_i A_i ) \\ D_i \sim LogNormal( \bar{D}_i, \sigma) \\ \bar{D}_i = \alpha + \sum_{k=1}^{K} \beta_k x_{i,k} \tag{5.2} \end{equation}\]
where \(A_i\) is the observed data measuring total settled area within a location \(i\). If the area is measured in hectares, then \(D_i\) represents people per hectare. \(\bar{D}_i\) is the expected population density on the log scale (i.e. the mean of the log-normal distribution), and \(\sigma\) is the residual variance term. \(K\) is the total number of covariates included in the model.
This hierarchical formulation has several advantages over the simple Poisson regression from Eq. (5.1):
- There is an added residual variance term \(\sigma\) that allows for over-dispersion of the Poisson,
- Covariates predict population density rather than counts, and
- The log-normal replaces the log-link function - acting as a stochastic log-link.
Over-dispersion means that the model is able to accommodate more residual variance in population counts than can be modelled with a Poisson distribution alone as given it’s lack of a variance parameter. In addition, by making the covariates predictors of population density rather than population counts the confounding effect of area is avoided. For example, two locations with identical covariate values and population densities could have very different population counts if the total amount of settled area is different. The hierarchical model explicitly accounts for this multi-level process.
Eq. (5.2) will serve as the core likelihood model for many of the model customizations described below.
5.2.5 Age-sex Structure
We can incorporate an age-structured sub-model (Boo et al. 2022a) if the household survey data contain counts \(M_{i,g}\) of people in each age-sex group \(g\) at each location \(i\). Such data allow for the estimation of population pyramids (i.e. proportions of the population in each age-sex group) and the production of age-sex-specific population estimates. A multinomial model can be added to Eq. (5.2) to achieve this:
\[\begin{equation} M_{i,g} \sim Multinomial(\theta_{r,g}, N_i) \\ \theta_{r,g} \sim Dirichlet(rep(1,g)) \tag{5.3} \end{equation}\]
where \(N_i\) is the total population at location \(i\) from Eq. (5.2). The population pyramid \(\theta_{r,g}\) is estimated independently for each region \(r\) with a flat Dirichlet prior. The Dirichlet prior enforces the assumption that individual elements of \(\theta_{r,g=1:G}\) are between zero and one and that they sum to one across all age-sex groups \(g\).
5.2.6 Random Intercept
Random effects are regression coefficients (e.g. \(\alpha\) and \(\beta\) above) that are dependent on other parameters. All of the regression coefficients shown above are fixed effects because they are not dependent on other parameters. Models that contain random effects are sometimes called mixed effects models because they contain fixed and random effects. Mixed effects models may have random intercepts, random slopes, or both.
An example of a random intercept in a population model is a a regression intercept \(\alpha\) (i.e. average population density) that is estimated separately for urban and rural areas in a way that accounts for the correlation between the two. Eq. (5.2) can be adjusted to include this random intercept \(\alpha_t\):
\[\begin{equation} N_i \sim Poisson( D_i A_i ) \\ D_i \sim LogNormal( \bar{D}_i, \sigma) \\ \bar{D}_i = \alpha_t + \sum_{k=1}^{K} \beta_k x_{i,k} \\ \alpha_t \sim Normal(\eta, \theta) \end{equation}\]
where \(t\) is the settlement type (i.e. urban or rural) that location \(i\) belongs to. \(\eta\) and \(\theta\) are the mean and standard deviation of \(\alpha\) among settlement types. The correlation between \(\alpha\) for the two settlement types is explicitly modelled because they are drawn from the same distribution, However these parameter estimates will still differ based on their fit to the data from each settlement type. This is a random intercept by settlement type, and it can help to account for the stratified sampling that household surveys often use to collect population data.
We can extend this concept to a random intercept by settlement type \(t\) and region \(r\) to account for additional spatial correlation where population densities from the same region are more similar to one another than population densities from different regions. Regions \(r\) can be defined as states or local government areas. This two-level random intercept \(\alpha_{t,r}\) (by settlement type and region) can be included as:
\[\begin{equation} N_i \sim Poisson( D_i A_i ) \\ D_i \sim LogNormal( \bar{D}_i, \sigma) \\ \bar{D}_i = \alpha_{t,r} + \sum_{k=1}^{K} \beta_k x_{i,k} \\ \alpha_{t,r} \sim Normal(\breve{\alpha}_{t}, \theta_{t}) \\ \breve{\alpha}_{t} \sim Normal(\bar{\alpha}, \eta) \tag{5.4} \end{equation}\]
where \(\breve{\alpha}_{t}\) and \(\theta_{t}\) are the mean and standard deviation (for each settlement type) of regression intercepts \(\alpha_{t,r}\) among regions. At the national level, \(\bar{\alpha}\) and \(\eta\) are the mean and standard deviation for \(\breve{\alpha}_{t}\).
This hierarchical random intercept can help to account for:
- Sampling that is stratified by settlement type, and
- Spatial autocorrelation within regions.
5.2.7 Hierarchical Variance
Similar to the hierarchical random intercept above, hierarchical variance by settlement type and region can also be used. This allows for uncertainty to be mapped and to see where residual variance is the greatest, providing more realistic ranges of uncertainty around population estimates in different regions and settlement types. Eq. (5.2) can be modified to have hierarchical variance \(\sigma_{t,r}\):
\[\begin{equation} N_i \sim Poisson( D_i A_i ) \\ D_i \sim LogNormal( \bar{D}_i, \sigma_{t,r}) \\ \bar{D}_i = \alpha + \sum_{k=1}^{K} \beta_k x_{i,k} \\ \sigma_{t,r} \sim HalfNormal( \breve{\sigma}_t, \theta_t ) \\ \breve{\sigma}_t \sim HalfNormal( \bar{\sigma}, \eta ) \tag{5.5} \end{equation}\]
Half-Cauchy distributions are also often recommended for modelling hierarchical variances rather than the Half-Normal that we have shown here (Gelman et al. 2013). Hierarchical variances can lead to convergence issues and care must be taken to specify priors that result in good convergence without being too influential on the posterior parameter estimates. It is often necessary to simplify the variance structure (e.g. fewer settlement types, or regions, or dropping one level entirely), especially if the sample size is low in some regions and/or settlement types.
5.2.8 Weighted-likelihood
Household surveys often implement a weighted sampling design known as PPS, or Probability Proportional to Size. This means that the probability of a location being selected for a survey is not random, it is dependent on the number of people (or households) within that area. Household surveys use weighted sampling to achieve a representative sample of households. If surveys were to use spatial random sampling, the results would likely bias towards rural areas because urban areas generally occupy less space on the landscape.
To use these datasets for population modelling, it is necessary to account for the bias that weighted sampling can introduce to avoid overestimating average population densities. Assuming sample weights \(w_i\) were used to collect a weighted sample of locations from a national sampling frame, a weighted-likelihood model can be developed that incorporates these weights to provide unbiased estimates of population densities. The first step is to calculate inverse weights and scale them to sum to one:
\[\begin{equation} m_i = \frac{w_i^{-1}}{\sum_{i=1}^I{w_i^{-1}}} \end{equation}\]
where \(I\) is the total number of observations used to fit the model. The scaled inverse weights \(m_i\) (or “model weights”) are used to weight individual samples in the likelihood by adjusting the variance term \(\sigma_i\):
\[\begin{equation} N_i \sim Poisson( D_i A_i ) \\ D_i \sim LogNormal( \bar{D}_i , \sigma_i ) \\ \sigma_i = \sqrt{ \frac{1}{m_i \theta^{-2}} } \tag{5.6} \end{equation}\]
where \(\theta\) is an estimated parameter that is a component of the variance, together with the model weights \(m_i\). The standard deviation for the log-normal \(\sigma_i\) in this model is location-specific, resulting in unbiased estimates of the mean and variance. This is because it gives more weight to the likelihood to locations that had lower probabilities of being included in the sample (e.g. for PPS household survey designs this would be locations with fewer people). The regression model for \(\bar{D}_i\) is not shown but it could be setup similar to Eq. (5.2).
The sample weights \(w_i\) are often unknown for unsampled areas where population predictions are needed. Because of this, a weighted average value for the variance term needs to be derived that is not location-specific:
\[ \bar{\sigma} = \frac{ \sum_{i=1}^I{ \sigma_i \sqrt{m_i} } } { \sum_{i=1}^I{ \sqrt{m_i} } } \]
This is a weighted average of \(\sigma_i\) across locations \(i\) - essentially factoring out the model weights \(m_i\). Model predictions of population density \(\hat{D}_i\) in locations where sampling weights \(w_i\) are unknown can be produced from:
\[ \hat{D}_i \sim LogNormal(\bar{D}_i, \bar{\sigma}) \]
5.2.9 Geostatistical Models
Geostatistics is a form of spatial statistics that explicitly model a continuous spatial phenomenon when observations are accurately georeferenced at particular sites (such as from a GPS location in a survey). Geostatistical models can help to estimate the outcome in unobserved locations, with the expectation that nearby locations are more similar than distant locations. While geostatistical modelling includes interpolation or smoothing methods such as Kriging, a model-based geostatistical approach (Diggle & Giorgi 2016a) makes it possible to incorporate spatial position into a statistical framework similar that shown in Eq. (5.1). The spatial information from the observations’ locations (in addition to observed covariate data) can improve the accuracy of population estimates. WorldPop has utilised geostatistical modelling across a number of applications - for example, producing population estimates for Cameroon, Papua New Guinea, Nigeria and the Democratic Republic of Congo (see chapter 16 and chapter 17), mapping the proportion of populations under 5 years of age (Alegana et al. 2015), producing high-resolution poverty estimates (Steele et al. 2017b), and estimating vaccination coverage (Utazi et al. 2019).
The general form of a model-based geostatistical framework is a mixed-effects regression model. This model includes fixed covariate effects plus a spatially correlated random effect for modelling spatial variation:
\[\begin{equation} N_i \sim Poisson( \mu_i ) \\ log(\mu_i) = \alpha + \beta x_i + Z(i) \tag{5.7} \end{equation}\]
where \(i\) indicates locations of observations, but these locations are taken as a spatial index within a fixed domain (\(s \in D \subset \mathbb{R}\)). \(Z(\cdot)\) is a spatially-continuous process that can be modelled as a Gaussian random field. Estimating the characteristics of a Gaussian random field is an important component of geostatistics, in particular the covariance (\(\Sigma\)) which describes how the dependence varies as a function of distance between the observations.
Geostatistical models are often implemented using the Integrated Nested Laplace Approximation (INLA) approach (Lindgren & Rue 2015, Krainski et al. 2018) because estimation of spatial covariance is very computationally-intensive for MCMC samplers. INLA is a more efficient alternative to traditional Bayesian MCMC sampling that provides rapid approximations of posterior probability distributions. Krainski et al. (2018) provide a good introduction to fitting Bayesian spatial models.
5.3 Conclusion
Bayesian statistical models used for bottom-up population estimates are powerful and versatile. This chapter summaries the core strategies that can be implemented to build appropriate statistical models for various data sets and applications. The methods outlined here are freely available to be built upon and customized for new applications.
Contribution
This chapter was written by Douglas Leasure, Andy Tatem and Chris Jochem