Integrating surveys with geospatial data through small area estimation to disaggregate SDG indicators: A practical application on SDG Indicator 2.3.1

1.Introduction

In an era characterized by the proliferation of new data sources and an unprecedented data revolution, the 2030 Agenda for Sustainable Development and the overall goal of leaving no-one behind (LNOB) have generated a tremendous increase in the demand of disaggregated data and statistics. In particular, in order to operationalize the overarching requirement of data disaggregation in the development of the Global SDG Indicator framework, the United Nations Statistical Commission postulated that “SDG Indicators should be disaggregated, where relevant, by income, sex, age, race, ethnicity, migratory status, disability and geographic location, or other characteristics, in accordance with the Fundamental Principles of Official Statistics”.

In this framework, traditional sample surveys implemented by National Statistical Offices (NSOs) can provide important information on the social, economic and environmental dimensions of target populations, representing the essential data source to produce the official estimates of about the 30% of Sustainable Development Goal (SDG) Indicators.11 However, these data sources alone are not enough to realize the ambitious goal of monitoring SDG Indicators by all relevant disaggregation dimensions and geographical areas. Indeed, despite collecting detailed information at relatively high frequency, most sample surveys are characterized by sample sizes that are either not large enough to guarantee reliable direct estimates for all sub-populations or that do not cover all possible disaggregation domains [1].

Issues of this kind can be addressed at different stages of the statistical production process. They can be tackled at the design stage, by adopting sampling strategies guaranteeing an observed set of sampling units for every disaggregation domain. Although potentially optimal, this approach normally results in an exponential increase of the sampling size and survey costs and complexity [2]. Furthermore, it is important to realize that, in practice, the anticipation of all possible future uses of survey data is virtually impossible, as “the client will always require more than is specified at the design stage” [3]. Alternatively, data disaggregation can be addressed at the data analysis stage, by adopting indirect estimation approaches borrowing strength from related disaggregation domains and/or time periods, thus resulting in an increase of the effective sample size [4]. Small area estimation (SAE) methods are among the possible indirect estimation approaches that can be adopted to deal with data disaggregation at the analysis stage. SAE techniques allow combining survey data with auxiliary information coming from additional data sources that are not affected by sampling error. Traditionally, SAE have relied on the integration of survey microdata with information from population and agricultural censuses or administrative records through explicit models linking the variable of interest to a set of auxiliary variables. However, with more and more data made available to National Statistical Systems (NSSs) from multiple innovative data sources, relying exclusively on auxiliary variables from traditional statistical sources for the production of small area estimates of SDG indicators is not considered as an efficient solution. In this respect, the 2030 Agenda explicitly stresses the need for new and enhanced data integration strategies, including the exploitation of the potential contribution to be made by geospatial information systems and other big data sources.

Within this framework, the present paper makes the case for the adoption of SAE and other indirect estimation methods to produce granular disaggregated estimates of SDG indicators, by integrating survey microdata with auxiliary information retrieved from “innovative” data sources, such as earth observation data. Indeed, relying on suitably implemented indirect estimation techniques such as SAE allows obtaining reliable disaggregated estimates of SDG indicators while managing survey costs and complexity. In particular, the integration of survey microdata with data from non-traditional sources offers the potential of producing timelier and more disaggregated statistics at higher frequencies than what allowed by traditional data sources alone. The paper is structured as follows. Distinguishing between area-level and unit-level SAE models, Section 2 presents a brief overview of these two approaches along with their relevant notation, context of usability, and potential source of auxiliary data. This section also discusses the main strengths and elements of cautions related to the use of geospatial auxiliary variables. Then, Section 3 presents some of the main fields of SAE application in official statistics, highlighting the few initiatives and examples of SAE for data disaggregation of official SDG indicators. To fill the gap of SDG relevant studies, Section 4 presents a practical case study on SDG indicator 2.3.1, measuring the average volume of production per labour unit of small-scale food producers, based on the Fay-Herriot (FH) [5] area-level model and combining the official integrated household survey of Mali with auxiliary variables retrieved from multiple geospatial information systems. The results of the case study highlight that, implementing the considered SAE approach, estimates precision is improved and predictions for out-of-sample domains can be produced. Finally, the main conclusions and ways forward are presented in Section 5.

2.Integrating survey data with additional data sources through small area estimation

Sample surveys, which are regarded as cost-effective means to collect detailed information at relatively high frequency over time, have a long history in the field of official statistics, and can be used to produce reliable estimates of parameters referred to total populations or to broad disaggregation domains. In this context, direct domain estimates of target parameters are statistics based solely on domain-specific sample data. Direct estimators are also known as design-based estimators, since they make use of sampling weights to produce inference on the target population [6]. One of the main requirements to achieve reliable disaggregated estimates by direct estimators is the presence of a sufficient domain sample size to yield adequate precision, or, in other words, a small estimated variance. When this circumstance is not verified, we are in the presence of so-called small-areas, i.e. disaggregation domains for which too little or no sampling observations are available [4]. It should be noted that, in practical statistical applications, it is quite rare to have an overall sampling size that is large enough to guarantee a sufficient number of observations for every possible disaggregation domain. Therefore, the use of indirect estimation techniques that “borrow strength” from auxiliary information on the population of interest [7] is often necessary. The range of possible approaches to produce indirect estimators is vast and goes from the implementation of design-based model-assisted approaches, such as the generalized regression estimator ([8, 9]) or the projection estimator ([10, 11]), to model-based approaches such as SAE ([4, 5, 12]). Contrarily to direct and model-assisted approaches, SAE model-based methods rely on explicit models and, consequently, the properties of resulting estimators are assessed under the adopted model assumptions. In particular, traditional SAE models are mixed models with area-specific random effects accounting for the variability between different areas not explained by auxiliary variables [4].

Although different in their specification, all SAE approaches share the same notation framework that is here introduced for the clarity of following sections. Let us consider a finite population U of N units that can be partitioned into D estimation domains U1,…,UD of sizes N1,…,ND. With d we denote the dth disaggregation domain, while i specifies the ith unit of the population. Let us now consider a random sample s∈S of size n (with S being the set of all possible sample s of size n that can be selected from U) and probability p⁢(s), the units of which can be used to produce direct estimates θ^d of target disaggregation parameters θd related to a variable of interest y. Typical examples of disaggregation parameters θd that are usually estimated for continuous variables are the domain total Yd=∑i∈Udyi and mean Y¯d=Yd/Nd. The well-known Horvitz-Thompson (HT) estimator of Yd and Y¯d can be expressed as Y^d=∑i∈sdwi⁢yi and Y¯^d=Y¯d/∑i∈sdwi, with wi=1/πi denoting the sampling weighs and πi=∑{s:i∈s}p⁢(s) the inclusion probability of unit i. It should be noted that Y¯^d has the functional form of a ratio estimator, as both its numerator and denominator are sampling estimates. The HT estimator of the total is design unbiased, while the one for the ratio is affected by a bias that tends to 0 with increasing values of nd. This means that their expected values are or tend to be equal to the parameter to be estimated [13]. As a consequence, their reliability is assessed only in terms of their precision, i.e. by the extent of their variance.

Direct HT estimators θ^d are usually characterized by unknown variance V⁢(θ^d) that needs to be estimated with adequate estimators v^⁢(θ^d), for a complete overview of which we refer to [9, 13, 14]. When the estimated variance is unacceptably high, SAE and other indirect estimation approaches can be used to increase estimates precision. Model-based SAE approaches allow considering the unexplained heterogeneity among domains, and have the potential of providing estimates that are more efficient than those produced with direct estimation methods. In addition, relying on SAE it is possible to predict indicator values also in out-of-sample domains. The literature on SAE classifies its models into two broad categories identified as area-level and unit-level models, which are briefly discussed in the two sections below. While area-level approaches relate a small area direct estimator θ^d to area-specific auxiliary information and can be adopted also when unit-level data is not available, unit-level models require access to microdata at the unit level, as they relate the unit values yi to unit-specific covariates [9].

Despite their increasing popularity, resorting to SAE should not be considered as the solution to any data disaggregation problem, and there are various considerations that NSOs should make before engaging in the production of indirect estimates. First of all, model-based approaches have stricter data requirements than direct estimation methods, with unit-level models being more data intensive than area-level ones. In this respect, the access to microdata on individual units may be limited by confidentiality concerns that need to be taken into account. Being based on models, after implementing SAE approaches, the underlying assumptions need to be carefully validated through adequate diagnostic techniques [27]. In addition, the bias of small area estimates needs to be measured to assess estimates reliability. This is generally done by means of the mean square error (MSE), which provides a combined indicator of estimates precision (variance) and accuracy (bias).

3.Use of small area estimation approaches for data disaggregation of SDG indicators

The empirical literature on SAE is very broad, with applications in many different fields of official statistics such as income and poverty, labour, health and agriculture. However, despite the great emphasis placed on data disaggregation in the context of the SDG monitoring framework and its overarching LNOB pledge,33 the number of examples of SAE techniques applied to official SDG indicators is still limited.

Being poverty mapping among the main applications of SAE, several case studies and references are available for the disaggregation of indicators related to Goal 1 on ending poverty. In particular, SAE techniques have been implemented to produce official sub-national estimates of SDG indicators 1.1.1 and 1.2.144 in countries such as Albania, Bolivia, Bulgaria, Cambodia, Chile, Ecuador, Indonesia, Mexico, Morocco and Sri Lanka ([23, 36]). Other applications relevant to income and poverty analysis, yet without a direct link with SDG indicators, can be found, for example, in Tanzania [37] and the United States [38].

Concerning Goal 2, aiming at ending hunger, achieving food security, improving nutrition and promoting sustainable agriculture, applications of SAE relevant to food security and malnutrition were found in Nepal [39], Ethiopia [40], and the United States [41]. However, the only application of indirect estimation techniques targeting specifically an indicator under Goal 2 was developed by the Food and Agriculture Organization of the United Nations (FAO) for indicator 2.1.2 on the prevalence of moderate and severe food insecurity in the population based on the food insecurity experience scale ([1, 11]). Concerning the agricultural component of Goal 2, evidence of empirical applications of SAE targeting SDG indicators under this goal were not found. Indeed, while SAE approaches have extensively been used to produce disaggregated estimates of crop yield and production measures (see [33, 34] for some examples), the use of indirect estimation approaches to produce disaggregated measures of agricultural labour productivity (indicator 2.3.1) or agricultural sustainability (2.4.1) are not a common practice.

Finally, a limited number of applications of SAE on indicators under Goal 4 [42], 5 [43], and 8 [44] were identified.

4.Empirical application of SAE on SDG Indicator 2.3.1 with the use of geospatial auxiliary variables

Target 2.3 of the 2030 Agenda for Sustainable Development aims to double the agricultural productivity and incomes of small-scale food producers by the end of the monitoring period. Progress towards the achievement of this target is monitored by two official SDG indicators, namely indicator 2.3.1 – measuring the average value of agricultural production per labour unit55 – and indicator 2.3.2 – estimating the average income from agricultural production activities of small-scale food producers. Indicator 2.3.1, which is the object of the presented case study, provides a measure of average partial factor productivity of agricultural holdings in a given year, and is currently disaggregated by the sex of the holding’s head and the size of the farm (small versus non-small). In particular, small-scale food producers are identified through an official definition developed by the FAO and endorsed by the Inter-Agency and Expert Group on SDG Indicators in September 2018 [45] in order to enhance international comparability. Although disaggregation at the subnational level is not among the mandatory disaggregation dimensions for reporting indicators under target 2.3, local estimates of indicators 2.3.1 and 2.3.2 may prove to be way more relevant than national aggregates for effective monitoring and decision making at the country level.

The typical data sources used to estimate SDG indicators on small-scale food producers are agricultural surveys, or household surveys integrated with modules on households’ agricultural activities. Being based on sample data, the production of reliable estimates of these two indicators at granular subnational level is usually not possible with standard design-based approaches, and –consequently – indirect estimation approaches need to be explored. As discussed in Section 3, the literature on small area income estimation is considerably wide, even if not necessarily targeting the estimation of incomes generated through agricultural production activities. Contrarily, the body of work on SAE approaches applied to indicators of labour productivity measures is still very little. In order to fill this gap, this section explores the application of a FH area-level model to produce small area estimates of SDG indicator 2.3.1 at the second administrative level (circles) of Mali, considering the integration of household survey data with area-level auxiliary information retrieved by multiple trustworthy geospatial information systems. Specifically, the presented SAE application is based on microdata from the Enquête Agricole de Conjoncture Intégrée aux Conditions de Vie de Ménages (EAC-I) 2017. The EAC-I is a multi-thematic cross-sectional household survey, implemented under the World Bank Living Standard Measurement Study (LSMS) programme, based on a nationally representative sample of about 8,390 households and with a specific focus on agriculture. In 2017, the sample units were divided into two groups, one of 3,813 households that received the full questionnaire, and one with remaining households that received a light version of the same questionnaire. For the purpose of this application, only the group of households that completed the full questionnaire could been considered, since this included the necessary variables to identify small-scale food producers and compute the targeted indicator. Considering that indicator 2.3.1 has a disaggregation dimension already embedded in its definition, i.e. the size of the farm, the sample that could be used to produce small area estimates for small scale food producers included only 1,637 households. Table 1 provides a summary of the sample size by region and circle, and the number of out-of-sample circles. In particular, the entire region of Kidal was left outside of the sample due to security reasons. In addition, the new region of Menaka had not been officially announced yet at the time of the survey and, for this reason, was not included in the sample.

Table 1 provides information on the sampling size by region and circles.

4.1Parameter of interest, selection of SAE approach, and considered geospatial auxiliary variables

Going back to notation introduced in Section 2, the average volume of production per labour unit to be estimated in each small area can be formalized as Y¯2.3.1,d=∑i=1Ndy2.3.1,iNd, with y2.3.1,i being the labour productivity of the ith small-scale food producer in circle d. The direct HT estimator of indicator 2.3.1 in the dth small area is Y¯^2.3.1,d=∑i∈sdwi⁢y2.3.1,i∑i∈sdwi, where the sampling weights are defined as in Section 2.

Table 2

Spatio-temporal resolution and sources of geospatial area-level covariates

Variable name	Spatial resolution	Temporal resolution	Source
Vol. fraction of coarse fragments	1 × 1 km	Static	ISRIC: World Soil Information
Nitrogen	1 × 1 km	Static
Sand	1 × 1 km	Static
Silt	1 × 1 km	Static
Clay	1 × 1 km	Static
Soil organic carbon	1 × 1 km	Static
Minimum temperature	4.5 × 4.5 km	Monthly	WorldCilm: Historical monthly weather data
Maximum temperature	4.5 × 4.5 km	Monthly
Precipitation	4.5 × 4.5 km	Monthly
Direct normal irradiation (Long-term yearly average)	0.3 × 0.3 km	1994-2018	Solargis
Diffuse horizontal irradiation (Long-term yearly average)	0.3 × 0.3 km	1994–2018
Air temperature (Long-term yearly average)	1 × 1 km	1999–2018
Vegetation indexes	5.5 × 5.5 km	Monthly	NASA EarthData
Elevation	1 × 1 km	Static	CGIAR CSI
Cropland	1 × 1 km	Annual	Zenodo
Bare ground	1 × 1 km	Annual
Built-up	1 × 1 km	Annual
Harvested area (major crops)	1 × 1 km	Annual	MAPSPAM
Production (major crops)	1 × 1 km	Annual

The accuracy of direct estimates, measured in terms of the estimated coefficient of variation (CV), was assessed against the same accuracy measure produced for EBLUP small area estimates obtained with the FH model in expression (1). This approach was selected in place of a unit-level method in order to produce a case study on SAE based on a model of simple implementation, only requiring access to area-level direct estimates and auxiliary information. In addition, as seen in Section 2.1, the EBLUP area-level estimator obtained from the model presented in expression (1) can be formulated as a linear combination of area-level direct and synthetic estimates, giving more weight to the former with increasing sampling size. Hence, using the FH SAE model can intuitively be seen as a way of improving direct estimates through a synthetic component based on external information. On the other hand, unit-level estimators, such as the one introduced in expression 2.2, do not take into account direct estimates and – as a consequence – the sampling design.

As area-level auxiliary variables for the implementation of the small area estimation model, various geospatial covariates were considered among the vast amount of publicly available candidates according to their potential capability of being good predictors for the average labour productivity in agriculture. In particular, covariates included in the first stage of selection were providing information on the following domains:

• • Soil characteristics: volume fraction of coarse fragments (> 2 mm), content of nitrogen, salt, silt, clary, and soil organic carbon.

• • Weather and climate: minimum and maximum temperature, precipitation quantity, direct normal irradiation, diffuse horizontal irradiation, air temperature, vegetation indexes.

• • Land cover: elevation, cover fraction of cropland, bare ground and extent of built up areas.

• • Harvested area and production of major crops (cotton, rice, sorghum, and wheat).

Table 2 presents the spatial and temporal resolution of each auxiliary variable along with the related source.

Table 3

Results of step-wise regression

Variable name	Unit of measure	lmg (%)
Production of cotton	(Metric ton)	24.0
Direct normal irradiation	(kWh/m2)	16.1
Production of wheat	(Metric ton)	14.6
Production of rice	(Metric ton)	11.7
Production of sorghum	(Metric ton)	11.4
Vol. fraction of coarse fragments	(%)	9.8
Soil organic carbon	(g kg-1)	8.7
Harvested area of rice	(hectare)	3.7

Values of considered geospatial predictor were initially available at the level of the cells of regular grids of different resolutions (spanning from 1 × 1 km to 5.5 × 5.5 km). Hence, being the basic FH approach based on auxiliary information referred to the small area of interest, data have been pre-processed in order to produce aggregates (totals or means depending on the variable) for the irregular polygons defining Mali’s circles.66

Figure 1.

Small area estimates of indicator 2.3.1 in Mali disaggregated by circle (second administrative division).

The initial set of potential predictors was then reduced adopting a stepwise regression, which was implemented using the area-level direct estimates of indicator 2.3.1 as dependent variable and the geospatial covariates as regressors.77 As result of the step-wise regression, only 8 auxiliary variables were retained (see Table 3) according to the Lindeman Merenda and Gold (LMG) factor, which represents a measure of the relative contribution of each predictor to the overall R square of the model. It is interesting to notice that most of the covariates considered as important by the selection approach provide information on either the quantity produced or the area harvested of Mali’s major crops. Other variables retained by the stepwise procedure measure the average direct normal irradiation, the average volume fraction of coarse fragments, and the average quantity of organic carbon in the soil.

4.2Smoothing variance of direct estimates

As seen in Section 2, among the necessary inputs to produce EBLUP small area estimates there are the estimated variances of direct estimates σ^ei2. Being based on few sampling observations, these estimates often need to be stabilized by means of some smoothing technique. For this case study, the approach based on the GVF with design effects [14] was adopted. Since the estimated design effects were far from being constant, the estimation domains were grouped into clusters with the objective of finding groups with similar design effects and inter-cluster correlations. The grouping was performed using the k-means clustering approach, which led to the identification of three groups with homogeneous design effects. The smoothed variances were used as input to the SAE model in place of the original variance of direct estimates. This process led to the elimination of the domain of Bamako, which was considered as out of sample due to an unacceptably high value of the smoothed variance.

Figure 2.

Boxplot of direct (left) and small area (right) estimates.

Figure 3.

Accuracy of direct and small area estimates and assessment of their linear relationship.

Figure 4.

Residuals and random effects of SAE model.

5.Conclusions and way forward

Monitoring the implementation of the 2030 Agenda for Sustainable Development and its overarching pledge to leave no one behind calls for more disaggregated data and SDG indicators than what available in most countries. In this context, sample surveys are the preferred data source for about the 30% of indicators in the SDG monitoring framework and can offer valuable information to measure the social, economic and environmental dimensions of sustainable development. However, traditional households and agricultural surveys are usually characterized by sampling sizes that are either too small to produce precise estimates, or that do not cover all disaggregation domains of interest. Hence, indirect estimation approaches such as SAE techniques can represent a valuable tool for NSOs and international organization to produce timely and granular disaggregated estimates of SDG indicators, allowing to contain the cost and complexity otherwise generated by the increase of sampling sizes. In particular, with the proliferation of new data sources such as geospatial and big data information systems, SAE models can be implemented by combining survey data with a vast amount of auxiliary information available at no or limited cost and at high frequency. In this respect, the body of literature and the number of case studies on SAE techniques applied to SDG Indicators can still be expanded. After a brief review of the main SAE approaches available along with their principal domains of application, this paper presents a case study based on the Fay-Herriot area-level SAE model to produce subnational estimates of SDG Indicator 2.3.1 on the average volume of production per labour unit obtained by small-scale food producers. This is done by integrating survey data with area-level auxiliary information retrieved from multiple geospatial information systems. The presented case study shows how the small area estimates of indicator 2.3.1 in Mali’s circles reach greater precision compared to direct estimates. In addition, adopting the considered indirect estimation approach, estimates for out of sample areas can also be produced.

The FH area-level model was selected in place of a unit-level method in order to provide a simple example of SAE based on an SDG indicator related to the agricultural sector development, only requiring access to area-level direct estimates and auxiliary information. In addition, using an indirect estimator – such as the area-level EBLUP – expressed as a linear combination of area-level direct and synthetic estimates, the SAE approach can intuitively be interpreted as a way of improving direct estimates through a synthetic component based on external information correlated with the phenomenon of interest. Future extensions of this study will compare the results obtained with the here considered area-level model with those produced by a unit-level approach. In this circumstance, both unit-level and sub-area (e.g. the enumeration area of the cell of a regular grid) level auxiliary variables will be considered as regressors.