Mapping poverty at multiple geographical scales

The first goal of the 2030 Agenda for Sustainable Development of the United Nations is to eradicate poverty in all its forms. Given that, current research trends are focusing on the complexity, spatial heterogeneity, and geographical roots of poverty to properly design out-of-poverty paths (Allard and Allard, 2017; Fan and Cho, 2021; Zhou and Liu, 2022). Poverty has a multifaceted nature compounded by its spatial dimension (Gauci, 2005): as an example, we can mention that widely used classifications divide it into urban and rural poverty (Christiaensen and Todo, 2014). Moreover, the study of spatial poverty traps, theorized as a persistent poverty status caused by location characteristics, e.g., remoteness, poor infrastructure and services, and/or excessive migration costs (Kraay and McKenzie, 2014) has received considerable attention.

For these reasons, poverty mapping is sparking interest in welfare economics and geography studies (Jean et al., 2016; Hall et al., 2023), as an essential methodology to support the investigation of the spatial distribution and regional characteristics of poverty. In addition, if performed at a granular level, it enables a geographical targeting, namely the area-driven allocation of resources for poverty alleviation (Elbers et al., 2007; Liu et al., 2017). This is recognized as an effective and cost-saving tool for poverty reduction if compared to individual targeting as the latter is characterized by high administrative costs of data collection and follow-up. Conversely, regional targeting does not require household tracking, easing the supervision and management processes (Bigman and Fofack, 2000; Galasso and Ravallion, 2005). Furthermore, it enables the diversification of policy instruments employed from area to area that can be easily combined with other antipoverty measures.

When performing poverty mapping, a natural problem relates to the choice of the spatial scale which strictly depends on the purpose of the analysis. This adds up to a more technical problem that arises when aggregating areal data from neighboring zones. Indeed, such an aggregation induces a smoothing effect that can mask or distort information on spatial heterogeneity (patchiness) of the target poverty measure (e.g., see Figure 1). This is known in geography as the scaling problem, and represents one of the aspects of the Modifiable Area Unit Problem (MAUP, Kolaczyk and Huang, 2001; Pratesi and Salvati, 2016). This effect strongly depends on the heterogeneity of the aggregated areas and might lead to misleading conclusions. A simple strategy to deal with it is to undertake the analysis by using multiple scales at once (Kolaczyk and Huang, 2001).

The choice of the spatial scale is also strongly connected with the availability of data sources. In particular, surveys on income, consumption and/or living conditions are commonly used in poverty mapping but typically provide reliable estimates at high levels of spatial aggregation (e.g., regional or national). Survey estimates are increasingly unreliable when considering finer spatial levels, translating into smaller and smaller sample sizes (Pratesi and Salvati, 2016). In addition, remote sensing (RS) or mobile phone usage data are recently employed in poverty mapping, especially in developing countries (Jean et al., 2016; Schmid et al., 2017). As opposed to survey data, such alternative sources are highly informative at fine spatial levels, while losing power when aggregated.

A good practice in poverty mapping is to promote the efficient integration of multiple sources by means of statistical models (Steele et al., 2017; Zhao et al., 2019). To do so, we have to optimize the informative power that distinct sources have at different scales while accounting for hierarchical dependencies among levels. In this sense, a multi-scale procedure enables such an integration, while overcoming the aforementioned scaling problem. To the best of our knowledge, no explicit formalization of a multi-scale modeling approach has been proposed in the poverty mapping literature. Some proposals exploit geo-statistical models to obtain a continuous surface estimation over the study area (Puurbalanta, 2021; Sohnesen et al., 2022), avoiding the scaling problem. Such methods, however, require the exact location of each household which is difficult to obtain due to data disclosure restrictions.

We propose a Bayesian multi-scale model for poverty mapping that combines survey and alternative sources, available at both finer and coarser spatial levels, and does not require the exact location of survey respondents. The quantity we target is a poverty rate, defined as the proportion of people whose income or living conditions indicator falls below a pre-specified threshold. Our methodology sets in the framework of the Small Area Estimation (SAE) literature, aiming at estimating population quantities in areas (or domains) not originally planned by the survey design (Tzavidis et al., 2018; Rao and Molina, 2015, Ch. 4), where area-specific sample sizes are likely to be too small to yield precise enough estimates, or even null. SAE methods have been widely used in poverty mapping (Molina et al., 2014; Casas-Cordero Valencia et al., 2016; Corral et al., 2022; Bikauskaite et al., 2022), and most of them rely on models defined at the unit level. Conversely, we adopt an area-level approach, namely integrating the area direct estimates with areal covariates. We remark that this approach requires neither survey microdata nor covariates defined for the entire population.

In this paper, poverty rates are modeled through a Beta-based model (Liu et al., 2014; Fabrizi et al., 2016b; Janicki, 2020). Specifically, we rely on the Extended Beta proposed in De Nicolò et al. (2022). Its features solve major issues encountered in poverty mapping such as the double-bounded support of poverty measures, the excess of zero/one values, and the intra-cluster correlation induced by the survey design. A new methodology to assess the uncertainty associated with out-of-sample areas is also proposed. We extend it to deal with multiple spatial scales, enabling us to overcome the scaling problem as well as to efficiently integrate multiple sources; in doing so, we provide a fresh contribution to poverty mapping literature.

Previous proposals that deal with multiple levels in SAE convey information in a single direction: from finer to coarser as the so-called sub-area or two-fold models (Torabi and Rao, 2014; Erciulescu et al., 2019), or from coarser to finer through geo-statistical models (Benedetti et al., 2022). Our proposal, instead, leverages the information streams in both directions by considering multiple predictors with shared random effects, introducing in the poverty mapping framework the approach proposed by Aregay et al. (2016, 2017) in the disease mapping literature.

A multi-scale approach must ensure the coherence of poverty estimates across levels: poverty rate estimates at a given low-resolution area must correspond to the population-weighted aggregation of higher-resolution area estimates. The misalignment among estimates obtained at multiple levels is settled by the so-called benchmarking techniques (Bell et al., 2013), reviewed in Section 4. Conventional benchmarking methods do not account for the bounded support of poverty rates and may return invalid outcomes, possibly lying below zero or above one. We propose a new benchmarking procedure that suitably restricts the parameter space to the unit interval and relies on the idea of posterior projection (Dunson and Neelon, 2003; Sen et al., 2018). This represents a contribution to the literature as it provides benchmarked estimates and related uncertainty measures that are fully consistent with the bounded support.

Our proposal of multi-scale modeling equipped with a suitable benchmarking procedure is applied to map poverty in Bangladesh at two different spatial layers: zilas, i.e. administrative districts, and upazilas, i.e. sub-districts comparable with counties or boroughs. We use data from the Bangladesh Demographic and Health Survey (DHS; Fabic et al., 2012) by considering as poverty measure the proportion of people in the area falling in the first quintile of the national distribution of the Wealth Index (WI; Pirani, 2014). WI is a composite score that summarizes the living conditions of a household and that can be read as a measure of socioeconomic status (Poirier et al., 2020). The insufficient survey sample sizes in most zilas and upazilas (one-third of the upazilas are out-of-sample) drive us to integrate survey with RS data. Open access RS data encompasses night light radiance, population structure and population density, along with geographical characteristics, land use, infrastructures, and other social and economic features known to be poverty predictors at granular levels (Zhou and Liu, 2022). Eventually, we compare the performance of our multi-scale model to alternatives by means of a design-based simulation study.

The paper is organized as follows. Section 2 introduces the notation; Section 3 is devoted to a detailed description of the proposed small area model defined at multiple levels, prior specification, and posterior inference algorithms, along with a review of other standard models. The benchmarking proposal is illustrated in Section 4. Section 5 presents a design-based simulation exercise that tests the performance of our model in comparison with alternatives. An application of poverty mapping in Bangladesh by using DHS and RS data is illustrated in Section 6. Conclusions are drawn in Section 7.

The statistical models discussed afterward primarily rely on survey data. In this section, we introduce the notation and basic survey quantities, such as the poverty rate estimators and related uncertainty measure, which constitutes the starting point of our modeling framework. Let $\mathcal{U}$ denotes a finite population of $N$ individuals, on which we define two nested partitions, although the following setting can be easily generalized to multiple partitions. The first partition divides $\mathcal{U}$ into $D$ non-overlapping domains $\mathcal{U}_{d}$ at a coarser level, labeled as areas, each constituted by $N_{d}$ units. A second partition at a finer level divides each $\mathcal{U}_{d}$, $d=1,\dots,D$, into $M_{d}$ sub-areas $\mathcal{U}_{dj}$ of size $N_{dj}$, so that $M=\sum_{d=1}^{D}M_{d}$ is the overall number of sub-areas.

A complex survey sample of overall size $n$ is drawn from $\mathcal{U}$ and, on its turn, it can be partitioned into area and sub-area sub-samples of sizes $n_{d}$ and $n_{dj}$, $\forall d,j$ respectively. The notations $\hat{\bar{Y}}_{d}$ and $\hat{\bar{Y}}_{dj}$ denote the estimators based on survey data (direct estimators) of a poverty rate for areas $d=1,\dots,D$ and the nested sub-areas $j=1,\dots,M_{d}$, respectively. Let $w_{dji}$ denotes the sampling weight of individual $i=1,\dots,n_{dj}$ pertaining to area $d$ and sub-area $j$, and $y_{dji}$ the indicator variable denoting the poverty status. As our target domains are not planned in the survey design, we employ an Hájek-type estimator (Hájek, 1971) defined as

	$\displaystyle\hat{\bar{Y}}_{dj}=\frac{\sum_{i=1}^{n_{dj}}w_{dji}y_{dji}}{\sum_{i=1}^{n_{dj}}w_{dji}},\quad j=1,\dots,M_{d},\ d=1,\dots,D$		(1)
	$\displaystyle\hat{\bar{Y}}_{d}=\frac{\sum_{j=1}^{M_{d}}\sum_{i=1}^{n_{dj}}w_{dji}y_{dji}}{\sum_{j=1}^{M_{d}}\sum_{i=1}^{n_{dj}}w_{dji}}=\frac{\sum_{j=1}^{M_{d}}\left(\sum_{i=1}^{n_{dj}}w_{dji}\right)\hat{\bar{Y}}_{dj}}{\sum_{j=1}^{M_{d}}\sum_{i=1}^{n_{dj}}w_{dji}}.$		(2)

Those estimators are approximately unbiased for the unknown population proportions denoted with $\theta_{d}$ and $\theta_{dj}$, assuming that $0<\theta_{d}<1$ and $0<\theta_{dj}<1$. Defining as $t$ the national poverty rate and as $s_{d}=N_{d}/N$ the known population share, it holds that $\sum_{d=1}^{D}s_{d}\theta_{d}=t$. Furthermore, the sub-area partition implies

$\displaystyle\sum_{j=1}^{M_{d}}s_{dj}\theta_{dj}=\theta_{d},\quad d=1,\dots,D,$

(3)

where $s_{dj}=N_{dj}/N_{d}$ with $\sum_{d}s_{d}=\sum_{j}s_{dj}=1$.

Survey estimates from (1) and (2) may be highly imprecise due to the small sample sizes. To assess the reliability of estimates, we need to introduce a measure of uncertainty, included also in the modeling framework of Section 3. To account for the complex design, we focus on the effective sample sizes $\tilde{n}_{d}$ and $\tilde{n}_{dj}$, i.e. the virtual sample sizes required by a simple random sample to retrieve estimators with standard errors equal to those of $\hat{\bar{Y}}_{d}$ and $\hat{\bar{Y}}_{dj}$, obtained under complex design. Focusing on the area case, we can express the sampling variance of $\hat{\bar{Y}}_{d}$ as $\theta_{d}(1-\theta_{d})\tilde{n}_{d}^{-1}$, being the estimator of a proportion. The equivalent sample sizes incorporate the gains and losses in efficiency due to stratification, weighting, and the possibly strong intra-cluster correlation. They can be expressed in terms of the design effect DEff_d, i.e. the ratio between the variance of the direct estimator and its variance under simple random sampling, as $\tilde{n}_{d}=n_{d}/\text{DEff}_{d}$. The same reasoning applies for $\hat{\bar{Y}}_{dj}$. We remark that, given the different correlation structures occurring within areas and sub-areas, generally $\tilde{n}_{d}\neq\sum_{j=1}^{M_{d}}\tilde{n}_{dj}$.

The estimated uncertainty measures allow the detection of spatial levels with unreliable estimates, thus requiring SAE techniques. In our framework, two different spatial levels are considered and we assume that they both present unreliable estimates as surveys are generally planned for national (or regional) aggregates.

Among SAE techniques, we focus on models defined at the area level. Small area models may be also defined at the unit (or individual) level and generally requires the individual auxiliary information to be known for the entire population. As novel data sources, such as RS data, provide areal information, the area-level models seem particularly suitable for poverty mapping.

We consider a hierarchical Bayesian model for proportions based on an extended Beta distribution, introduced in De Nicolò et al. (2022), by extending it for multiple spatial scales. We recall its single-level specification in Section 3.1, whereas multi-level specifications follow in Section 3.2.

The definition of area level estimates through (3) preserves the coherency of poverty rates between the two levels. However, the coherency with respect to higher levels of aggregation (e.g., national) is not guaranteed. Indeed, the poverty rate at higher levels may be known in population or reliably estimated through surveys; such value, hereafter called benchmark, might be used to constrain the linear combination of sub-area estimates. In this section, we consider the problem of finding a set of constrained estimators $\tilde{\bm{\theta}}=(\tilde{\theta}_{11},\dots\tilde{\theta}_{DM_{d}})^{T}$, for the target parameters $\bm{\theta}=(\theta_{11},\dots\theta_{DM_{d}})^{T}$. The linear constraint is imposed by the poverty rate at the finest spatial scales and it is defined by the combination of both equations in (3) as

$\displaystyle\sum_{d=1}^{D}\sum_{j=1}^{M_{d}}q_{dj}\tilde{\theta}_{dj}$

$\displaystyle=$

$\displaystyle t,$

(10)

where $q_{dj}=s_{dj}s_{d}$ denotes the population share of sub-area $dj$ at national level. Note that $\sum_{d}\sum_{j}q_{dj}=1$, and $t$ is, in our case, the poverty rate at the national level.

The problem of benchmarking has been tackled from different perspectives in the Bayesian literature. The most popular strategy follows a decision-theoretic framework: starting from a loss function $L(\bm{\theta},\tilde{\bm{\theta}})$, the benchmarked estimators are obtained by minimizing the posterior risk $\mathbb{E}[L(\bm{\theta},\tilde{\bm{\theta}})|\mathbf{y}]$ under linear constraint (Datta et al., 2011). Despite its simple implementation, this method is applied to point estimators obtaining only a set of benchmarked point estimators. Thus, it represents a non-fully Bayesian strategy, not delivering posterior measures of uncertainty. To solve this problem, fully Bayesian methods have been developed. Among the others, Zhang and Bryant (2020) include the constraint through a suitable prior distribution; Janicki and Vesper (2017) search for a constrained joint posterior distribution with minimum Kullback-Leibler distance to the unconstrained one, whereas Okonek and Wakefield (2022) opt for an importance sampling approach.

Our approach relies on the concept of posterior projection (Dunson and Neelon, 2003; Sen et al., 2018) that was already considered by Patra (2019) within the SAE framework. In this case, the posterior samples from $\theta_{dj}|\mathbf{y}$ are projected into the feasible set defined by the benchmarking constraint, minimizing the distance from the original ones. Specifically, this distance is defined by means of a loss function. In contrast to other fully Bayesian methods, such an approach is more computationally efficient, in particular when the fitted model is not trivial as in our case. The choice of a suitable loss function for our inferential problem is discussed in Section 4.1 and the main results are illustrated in Section 4.2.

In this section, we study how statistical models of Section 3 deal with a multi-scale data problem by evaluating the frequentist properties of resulting estimators. To this aim, we reproduce a framework in which the goal is making inference on a proportion and, hence, a dichotomous response is generated for the whole synthetic population. We opt for a design-based simulation study, being able to reproduce the peculiarities of a multi-scale problem within a survey design setting. To this aim, we consider two aggregation levels, i.e., areas and sub-areas, and two different scenarios that contemplate diverse ratios of between and within area variability.

For each scenario, a population of $N=180,000$ individuals grouped into $3,600$ clusters is generated, mimicking the two-stages structure characterizing the scheme of several common surveys. Units are divided into $D=30$ areas and $M=150$ sub-areas, i.e., each area has $M_{d}=5$ sub-areas. The whole procedure for generating the population and sample drawing is fully detailed in Section S3 of the Supplementary Material. A continuous variable is generated from a three-fold Gaussian model with cluster effect scale $\sigma_{c}=0.2$, area effect scale $\sigma_{a}$ and sub-area effect scale $\sigma_{s}$. In scenario 1, the sub-area variation ($\sigma_{s}=0.13$) prevails on the area variation ($\sigma_{a}=0.08$), inducing a marked heterogeneity within the areas. In scenario 2, the area variation ($\sigma_{a}=0.13$) prevails on sub-area variation ($\sigma_{s}=0.08$), leading to homogeneous areas. Such values for the scale parameters are retrieved from the variance decomposition related to real data discussed in Section 6. Lastly, the dichotomous response is obtained by assigning value 1 to individuals below the first quintile of the generated continuous variable, and 0 otherwise.

From each synthetic population, $B=1,000$ samples are drawn following a two-stages sampling procedure with rates of approximately 2%. To account for the presence of out-of-sample sub-areas, each drawn sample does not include units from 37 sub-areas belonging to 15 areas. In this way, half of the areas include individuals from all the sub-areas, and, in the other half, units pertaining to two or three sub-areas are not observed.

The frequentist properties of the benchmarked estimators $\hat{\tilde{\theta}}_{dj}$ and $\hat{\tilde{\theta}}_{d}$ obtained under different models are compared: S-MS, SA and I-MS. Note that, with the suffix -Aggr we denote area-level estimates defined as in (3). Let us indicate with $\theta_{i}$ as the true value for the generic area or sub-area $i$ and $\hat{\tilde{\theta}}_{i}^{(b)}$ as the benchmarked estimate at iteration $b$. To determine the 90% credible intervals we consider 5^th and $95^{th}$ posterior percentiles, labeled as $\tilde{\theta}_{i,L}^{(b)}$ and $\tilde{\theta}_{i,U}^{(b)}$. We define bias, root mean squared error (RMSE) and coverage for 90% credible intervals for each area or sub-area as

		$\displaystyle\text{Bias}\left[\hat{\tilde{\theta}}_{i}\right]=\frac{1}{B}\sum_{b=1}^{B}\left(\hat{\tilde{\theta}}_{i}^{(b)}-\theta_{i}\right),\quad\text{RMSE}\left[\hat{\tilde{\theta}}_{i}\right]=\sqrt{\frac{1}{B}\sum_{b=1}^{B}\left(\hat{\tilde{\theta}}_{i}^{(b)}-\theta_{i}\right)^{2}},$
		$\displaystyle\text{Cov}\left[\tilde{\theta}_{i}\right]=\frac{1}{B}\sum_{b=1}^{B}\bm{1}\left\{\theta_{i}\in\left[\tilde{\theta}_{i,L}^{(b)};\tilde{\theta}_{i,U}^{(b)}\right]\right\}.$

As summary measures of bias and RMSE, we define the mean absolute percentage mean squared error (MAPE) and the average relative root mean squared error (ARRMSE) as

$$\text{MAPE (\%)}=\frac{100}{BI}\sum_{b=1}^{B}\sum_{i=1}^{I}\frac{\hat{\tilde{\theta}}_{i}^{(b)}-\theta_{i}}{\theta_{i}},\quad\text{ARRMSE}=\frac{1}{I}\sum_{i=1}^{I}\frac{\text{RMSE}\left[\hat{\tilde{\theta}}_{i}\right]}{\theta_{i}}.$$

Lastly, as a measure of the model goodness-of-fit, the leave-one-out information criterion (LOOIC, Vehtari et al., 2017) is considered.

Interesting cues about the performances of the compared area-specific estimators can be deduced from Figure 2. It reports the relative RMSEs, obtained by dividing each RMSE by the corresponding S-MS-Aggr RMSE, jointly with biases and frequentist coverages versus the population values. In both scenarios, the S-MS-Aggr estimators are characterized by lower biases, if compared to I-MS and SA cases: such differences become relevant for areas whose true values are far from their average (about 50% of them). Such a decrease in the bias leads to some differences also in the behaviour of RMSE too: S-MS-Aggr produces more efficient estimators in the tails of the ensamble distribution of population poverty rates $\{\theta_{1},\dots,\theta_{D}\}$, becoming less efficient in the central part of the distribution. By focusing on the 90% credible intervals, we observe that S-MS-Aggr is affected by a slight under-coverage (median: 0.83 in scenario 1 and 0.81 in scenario 2). However, the performances deteriorate more slowly if compared to the other methods in the non-central parts of the ensamble distribution. Table 1 summarises the overall performances of the methods across areas. The MAPE points out the ability of S-MS-Aggr in reducing the bias, registering lower values both in areas with and without out-of-sample sub-areas; conversely, ARRMSE highlights that the average performances are quite similar.

We can point out that the I-MS model, characterized by the absence of an area-specific random effect, leads to poorer performances if compared to all the other strategies. This is more evident in Scenario 2 where the data-generating process mimics the case $\sigma_{a}>\sigma_{s}$. Targeting the estimation of sub-area parameters, summaries reported in Table S1 of the Supplementary Material point out that the three models behave quite similarly, with the exception of I-MS, which still shows lower performances in Scenario 2. Lastly, we highlight that the average LOOIC of the S-MS model is widely lower both for sub-area (Table S1) and area (Table 1) levels in both scenarios, denoting a remarkable improvement in the predictive abilities of this model.

In summary, the estimators under the S-MS model have more stable performances along the support. They are characterized by a more moderate shrinkage and, for this reason, they perform better on the tails of the poverty rates ensemble distribution, and slightly better than their competitors on average. Small and, especially, large poverty rates are relevant for the implementation of geographically targeted alleviation policies as they can identify poverty hotspots. Therefore, an improvement in the estimation of such values may be of great importance in this context.

The aim of this section is to map poverty in Bangladesh at multiple levels of administrative divisions: the Administrative Level AL-2, comprising $D=64$ zilas as areas, and the AL-3 with $M=544$ upazilas as sub-areas. The target poverty measure is the proportion of people below the 20th percentile of the WI national distribution. We consider DHS survey data complemented with RS and geographical data available from multiple open sources. Data and sources are shortly illustrated in Section 6.1, whereas results are discussed in Section 6.2.

In this paper, we introduce a multi-scale modeling framework for poverty mapping at multiple spatial resolutions to avoid the scaling problem. The main aim is to produce reliable poverty maps that are coherent at different geographical layers. Such tools are valuable for poverty evaluation, poverty targeting and to develop place-based policies. On this line, we complement the proposal with a novel benchmarking algorithm that ensures the accordance of poverty estimates across levels. This allows to preserve poverty rate estimates and associated credible intervals within their support, restricted to $(0,1)$.

The simulation study we carried out highlights the effectiveness of our proposal, showing better performances with respect to existing alternatives. In particular, this is evident for poverty rates laying in the tails of their distribution, which are less shrunken towards the average. In this way, our modeling strategy is able to retain domains with extremely high poverty rates: this has relevant practical implications since poverty hotspots need to be identified as being the target of poverty relief actions. In our application on Bangladeshi data, rural and remote regions that mostly overlap with droughts- and floods-prone zones with high exposure to climate change effects, are classified as very poor. In this sense, geo-targeted policies may foster risk mitigation strategies (e.g., insurance) as well as further expand irrigation coverages via specific water management policies (Pal et al., 2011; Hossain et al., 2022). Furthermore, our multi-scale poverty mapping approach permits us to spot critical zones located within larger administrative units.

It is possible to extend our proposal to explicitly account for the spatial structure in the statistical model, similarly to Aregay et al. (2016). However, results of our application do not provide evidence of a residual spatial trend. This can be due to the spatial informative power of remote sensing covariates employed in the model. Another future extension consists in considering also other quintiles of the WI index national distribution through a multivariate framework. Moreover, moving beyond poverty measurement, our proposal may be applied to map other well-being measures such as inequality indicators.

Work supported by the Data and Evidence to End Extreme Poverty (DEEP) research programme. DEEP is a consortium of the Universities of Cornell, Copenhagen, and Southampton led by Oxford Policy Management, in partnership with the World Bank - Development Data Group and funded by the UK Foreign, Commonwealth & Development Office. The work of Silvia De Nicolò was partially funded by the ALMA IDEA 2022 under Grant J45F21002000001 supported by the European Union - NextGenerationEU. The work of Aldo Gardini was partially supported by MUR on funds FSE REACT EU - PON R&I 2014-2020 and PNR (D.M. 737/2021) for the RTDA_GREEN project under Grant J41B21012140007.