Spatio-Temporal small area estimation of subdistrict per capita expenditure: Integrating panel survey and geospatial big data

Abstract

The availability of accurate poverty data at the granular (subdistrict) level remains a significant challenge for official statistics in developing countries, due to limited survey samples, creating information “blind spots” for policy-making. This study aims to address this issue by estimating subdistrict per capita expenditure in Medan Municipality and Langkat Regency by integrating official survey data with geospatial big data. Utilizing a Spatio-Temporal Hierarchical Bayesian Small Area Estimation (ST HB SAE) framework, this research leverages the 2021–2024 National Socioeconomic Survey (Susenas) panel data structure and auxiliary variables from satellite imagery, including Night Time Light (NTL), vegetation index (NDVI), built-up index (NDBI), and air pollutants (NO2). The developed ST HB SAE model incorporates spatial effects (Leroux CAR) and temporal autoregressive (AR1) processes to enhance estimation precision. Results demonstrate that this model significantly outperforms direct estimation and nonspatio-temporal models in terms of Root Mean Squared Error (RMSE) reduction and Coefficient of Variation (CV) stability. Geospatial variables are strongly correlated with welfare indicators, enriching information in sample-deficient areas. Furthermore, a benchmarking process ensures that the estimates are consistent with official regency/municipality level aggregates and meet official statistical standards for operational adoption. This approach offers a cost-effective, robust solution for National Statistical Offices (NSOs) to monitor welfare dynamics in small areas, supporting evidence-based development planning.

Keywords

Official statistics small area estimation Spatio-Temporal hierarchical Bayesian geospatial big data

1 Background

Poverty remains a primary global concern.¹ SDG-1 affirms the international commitment to end poverty in all its forms by formulating targeted policies based on accurate data.² However, National Statistical Offices (NSOs) in many developing countries face a major challenge: the high demand for micro-level data often clashes with limited survey budgets. This measurement challenge in granular domains renders policies ineffective and biased.³

In Indonesia, although the national poverty rate declined to 9.03 percent in March 2024,⁴ regional welfare inequality remains palpable, particularly between urban and rural areas.⁵ This inequality reflects structural issues influenced by economic, spatial, social, and environmental factors.^6,7 Therefore, poverty measurement that accounts for spatial and temporal variations becomes crucial to observe the dynamics of inequality and policy effectiveness across periods.

Statistics Indonesia (BPS) measures poverty based on average per capita expenditure as a proxy for household economic welfare. Measurement is conducted by comparing per capita expenditure against the poverty line, which reflects the minimum basic needs of households to achieve a decent standard of living.⁸ Based on this comparison, households with expenditures falling below the poverty line are categorized as poor.

The National socioeconomic Survey (Susenas) serves as the primary source of poverty data in Indonesia.⁹ However, due to sample size limitations, direct estimation from Susenas at the subdistrict level tends to yield high variance and low reliability. Consequently, official poverty data at the subdistrict level remain unpublished, creating an information “blind spot” for local governments that require precise data for allocating social assistance programs. This condition has driven various statistical studies to seek alternative approaches that increase estimation precision without increasing survey burden.

Small Area Estimation (SAE) is a widely used method to produce reliable estimates for more granular areas. This method improves estimation precision by combining information from other variables through a “borrowing strength” mechanism.^10,11 Various studies have demonstrated the effectiveness of SAE in reducing poverty estimation variance in small areas without increasing survey costs.¹² Nevertheless, most SAE models currently used assume that each location is independent and static. In reality, welfare phenomena exhibit interregional dependence and intertemporal connections.

SAE studies in Indonesia generally rely on administrative data such as Village Potential (PODES) or census data as auxiliary variables.¹³ Although easily accessible, these data are not available annually. PODES data are collected three times within ten years, while census data are updated every ten years. Such conditions render the data unable to reflect current field conditions. The limitations of census and administrative data hamper the model's ability to detect socioeconomic dynamics in the short term.¹⁴ Therefore, alternative auxiliary variables that are more dynamic and can cover small areas continuously are needed.

To fill this methodological gap, this study proposes an approach to modernize official statistics by integrating the Susenas Panel and Geospatial Big Data. Advances in remote sensing technology and geospatial big data processing offer a significant opportunity to enrich SAE models. Geospatial big data indicators have proven to correlate strongly with economic activity and community welfare.^15–30 For instance, Night Time Light (NTL) data used in Uganda were able to predict poverty,²⁶ while Normalized Difference Vegetation Index (NDVI), Normalized Difference Built-up Index (NDBI), and Modified Normalized Difference Water Index (MNDWI) used in Guangzhou were utilized to predict per capita expenditure.³⁰ Meanwhile, in Indonesia, NTL, NDVI, NDWI, Land Surface Temperature (LST), Carbon Monoxide (CO), Nitrogen Dioxide (NO2), and Point of Interest (POI) data can be used to form multidimensional poverty indices.¹⁶

Such data can be obtained quickly, are relatively inexpensive, and can produce sufficient information to represent the population,³¹ making them ideal as auxiliary variables in small area welfare estimation. Furthermore, the use of the Susenas Panel design from 2021 to 2024 offers a unique advantage, as data are collected from the same group of respondents repeatedly. By leveraging this dimension, SAE models can be developed into Spatio-Temporal SAE, where spatial and temporal variations are analyzed simultaneously to produce estimates that are more stable, adaptive, and policy-oriented.

Although the potential for integrating official survey data with geospatial big data has been explored, its application in the Indonesian context remains limited. Most previous studies focused on EBLUP SAE models without accounting for spatial and temporal structures,^15,19 composite index creation,^16,18 or directly predicting poverty.¹⁷ Few studies have applied the Spatio-Temporal Hierarchical Bayesian SAE (ST HB SAE) framework, whereas this approach allows for more realistic estimation by accommodating inter-regional and inter-period dependence. This limitation indicates a methodological and empirical gap that can be filled by the ST HB SAE approach, especially with the support of Susenas panel data and geospatial big data based auxiliary variables.

This study aims to estimate the average per capita expenditure in every subdistrict in Medan Municipality and Langkat Regency using the ST HB SAE model, which combines the strength of official panel survey data and geospatial big data variables. Methodologically, this study develops the implementation of ST HB SAE based on hierarchical modeling to produce spatially and temporally consistent posterior estimates of average per capita expenditure. Additionally, a benchmarking process is conducted to ensure the consistency of model estimation results with official aggregate figures at the regency/municipality level. This step is crucial so that estimation results can be operationally adopted into the BPS official statistics system. Practically, the research results are expected to support the modernization of BPS official statistics and strengthen evidence-based development planning at the subdistrict level.

2 Data used

This study integrates official survey data with geospatial big data indicators. The scope of the study covers 44 subdistricts in Medan Municipality (representing urban areas) and Langkat Regency (representing rural areas) over the period of 2021–2024. Medan Municipality and Langkat Regency were selected to represent the distinct urban-rural dichotomy of North Sumatra, as both regions face persistent welfare challenges and slow poverty reduction despite their contrasting economic structures. This integration of survey data with geospatial big data aims to enhance small area estimation while maintaining consistency with the official statistics framework of BPS.

2.1 Official survey data (Susenas panel 2021–2024)

The primary data source utilized in this study is Susenas, conducted by BPS. Unlike conventional cross-sectional approaches, this research leverages the structure of the 2021–2024 Susenas Panel data. The use of panel data enables the model to capture the dynamics of welfare changes within the same households over time, providing crucial temporal information for poverty estimation.

The target variable is the average per capita expenditure (on a logarithmic scale), calculated as total household consumption expenditure on food and nonfood items divided by the number of household members. Direct estimates for each subdistrict are calculated using the following formula:

\begin{aligned} {\hat{Y}}_{i t}^{d i r} = \frac{\sum_{j = 1}^{n_{i t}} w_{i j t} y_{i j t}}{\sum_{j = 1}^{n_{i t}} w_{i j t}} \end{aligned}

(1)

where

{\hat{Y}}_{i t}^{d i r}

presents the direct estimate of the average per capita expenditure in subdistrict i at year t;

w_{i j t}

is the sampling weight of the j -th household in subdistrict i at year t according to the Susenas design; and

y_{i j t}

is the per capita expenditure of the j -th household in subdistrict i at year . t The weights used to calculate the average per capita expenditure are based on population projection results from the 2015 Inter Census Population Survey 2015 (SUPAS 2015) (BPS, 2023a).

2.2 Geospatial big data

This study extracts various geospatial big data indicators to serve as auxiliary variables. The selection of variables is grounded in empirical studies demonstrating strong correlations between geospatial big data indicators and economic activity and social welfare.^{11,15,17,26,30}

All datasets were accessed via the Google Earth Engine (GEE) and OpenStreetMap platforms. Each variable was extracted as medium resolution rasters, followed by pre-processing steps including cloud masking, median compositing, and radiometric normalization to ensure quality and interannual consistency (Table 1). To maintain temporal alignment with the observations, all imagery was compiled during the March–July period of each year, coinciding with the Susenas data collection schedule. Subsequently, each raster layer was aggregated to the subdistrict level by calculating the mean value within administrative boundaries using official BPS shapefiles and the WGS-84 (EPSG 4326) projection system. The processing yielded a spatio-temporal panel dataset comprising six geospatial variables (Table 1) across 44 subdistricts (21 in Medan Municipality and 23 in Langkat Regency) over four observation years (2021–2024).

Table 1.
Description of geospatial big data sources, variables, and preprocessing specifications.

Source Preprocessing Variable Band Unit

NOAA-VIIRS Bidirectional Reflectance Distribution Function (BRDF) Correction NTL nanowatts/cm²/sr

Sentinel-2 Harmonized Low cloud image selection, cloud and shadow masking, temporal gap filling. NDVI Band 8 (NIR) dan Band 4 (Red) Index

NDBI Band 11 (SWIR) dan Band 8 (NIR) Index

MNDWI Band 3 (Green) dan Band 11 (SWIR) Index

Modis Conversion from Kelvin to Celsius; daily compositing and annual averaging. LST °C

Sentinel-5p Masking of zero/negative values; annual temporal compositing. CO Mol/m²

NO₂ Mol/m²

Open Street Map POI extraction POI Density amenity, healthcare, shop Count

Calculation of average Euclidean distance from grid center to nearest facility. Distance to nearest POI amenity, healthcare, shop Meters

Source	Preprocessing	Variable	Band	Unit
NOAA-VIIRS	Bidirectional Reflectance Distribution Function (BRDF) Correction	NTL		nanowatts/cm²/sr
Sentinel-2 Harmonized	Low cloud image selection, cloud and shadow masking, temporal gap filling.	NDVI	Band 8 (NIR) dan Band 4 (Red)	Index
NDBI	Band 11 (SWIR) dan Band 8 (NIR)	Index
MNDWI	Band 3 (Green) dan Band 11 (SWIR)	Index
Modis	Conversion from Kelvin to Celsius; daily compositing and annual averaging.	LST		°C
Sentinel-5p	Masking of zero/negative values; annual temporal compositing.	CO		Mol/m²
NO₂		Mol/m²
Open Street Map	POI extraction	POI Density	amenity, healthcare, shop	Count
Calculation of average Euclidean distance from grid center to nearest facility.	Distance to nearest POI	amenity, healthcare, shop	Meters

3 Methodology

The research methodology is designed through five systematic stages: (1) selection of auxiliary variables, (2) variance smoothing, (3) Spatio-Temporal Hierarchical Bayesian modeling, (4) parameter estimation and inference, and (5) evaluation and benchmarking.

3.1 Selection of auxiliary variables

The initial step is to select geospatial big data variables that exhibit a strong linear relationship with per capita expenditure. The selection was conducted using Pearson correlation ( $r$ ). The Pearson correlation coefficient between two variables is formulated as follows³²:

\begin{aligned} r = \frac{\sum_{i = 1}^{m} \sum_{t = 1}^{T} (x_{i t} - \bar{x}) (y_{i t} - \bar{y})}{\sqrt{\sum_{i = 1}^{m} \sum_{t = 1}^{T} {(x_{i t} - \bar{x})}^{2}} \sqrt{\sum_{i = 1}^{m} \sum_{t = 1}^{T} {(y_{i t} - \bar{y})}^{2}}} \end{aligned}

(3)

Where r denotes the Pearson correlation coefficient, $x_{i t}$ represents the auxiliary variable, $y_{i t}$ represents the log direct estimate of per capita expenditure, m denotes the number of subdistricts and T the number of years. In this study, only variables exhibiting moderate to strong correlation (|r|≥0.4)) were retained as predictors in the model.^33,34 This moderate to strong correlation threshold was chosen as an initial screening criterion to ensure a reasonable signal to noise ratio in the hierarchical bayesian framework. According to Evans,³⁵ a correlation coefficient in this range represents a meaningful relationship in socioeconomic studies. Importantly, within a hierarchical bayesian model, the final contribution of auxiliary variables is further regularised through coefficient shrinkage. This approach allows the model to borrow strength from relevant geospatial indicators without introducing excessive complexity or overfitting that might arise from including weaker predictors.

3.2 Variance smoothing

Direct estimates from Susenas at the subdistrict level often exhibit unstable variances due to small sample sizes. To address this issue, the Generalized Variance Function (GVF) approach is employed, with the form³⁶:

\begin{aligned} V a r ({\hat{Y}}_{i t}^{d i r} | θ_{i t}) = \frac{f (θ_{i t})}{n_{i t}} \end{aligned}

(4)

Where $n_{i t}$ denotes the sample size and $f (θ_{i t})$ represents the variance function. The specified variance function is defined as:

f (θ_{i t}) = θ_{i t}^{2}

Thus, the smoothing model is written in equation (5) as follows:

\begin{aligned} {\hat{V}}_{i t}^{d i r} = ψ \frac{{({\hat{Y}}_{i t}^{d i r})}^{2}}{n_{i t}} + ε_{i t} \end{aligned}

(5)

With $ψ = 1 / d e f f$ . Parameter estimation is performed using Generalized Least Squares (GLS). The GLS approach is preferred because survey data with small domains tend to exhibit varying variances across areas.

3.3 Spatio-Temporal HB model specification

This study extends the Fay-Herriot (FH) model into a Spatio-Temporal Hierarchical Bayesian framework by incorporating the likelihood, the linking model, and spatial and temporal random effects. The likelihood is defined as in equation (6):

\begin{aligned} y_{i t} | θ_{i t}, {\tilde{V}}_{i t}^{(log)} \sim N (θ_{i t}, {\tilde{V}}_{i t}^{(log)}) \end{aligned}

(6)

Where $y_{i t} = \log ({\hat{Y}}_{i t}^{d i r})$ , ${\hat{Y}}_{i t}^{d i r}$ representing the log of the direct estimate for subdistrict i at time t; $V_{i t}$ represents the result of variance smoothing and transformed to the log scale via the delta method. The value of $μ_{i t}$ is linked to big data based auxiliary variables, year effects, and spatio-temporal effects through the following linking model:

\begin{aligned} θ_{i t} = x_{i t}^{T} β + δ_{t} + s_{i t} \end{aligned}

(7)

Where $x_{i t}$ denotes the vector of auxiliary variables used, $β$ represents the fixed coefficients measuring the influence of the auxiliary variables, $δ_{t}$ indicates the specific year effects, and $s_{i t}$ is the random effect encompassing spatial and temporal variations. To capture the spatial dependence between adjacent subdistricts, the Proper Conditional Autoregressive (CAR) model proposed by Leroux is employed.³⁷ The spatial effect is expressed in equation (8) as follows:

\begin{aligned} z_{t} \sim C A R_{p r o p e r} (τ_{s}, ρ_{s p}), t = 1, \dots, T, \end{aligned}

(8)

The proper CAR distribution is defined with respect to the adjacency matrix $W = (w_{i j})$ , where $w_{i j}$ = 1 subdistricts i and j share a boundary and $w_{i j} = 0$ otherwise and $n_{i} = \sum_{j \neq i} w_{i j}$ denotes the number of neighbours of area i. The corresponding full conditional distribution is:

\begin{aligned} z_{i t} | z_{- i, t}, τ_{s}, ρ_{s p} \\ \sim N (\frac{ρ_{s p} \sum_{j \neq i} w_{i j} z_{j t}}{ρ_{s p} n_{i} + (1 - ρ_{s p})}, \frac{1}{τ_{s} [ρ_{s p} n_{i} + (1 - ρ_{s p})]} \end{aligned}

(9)

Where $ρ_{s p} \sim U (0.05, 0.95)$ , $σ_{s} \sim H a l f - N o r m a l (0, {2.5}^{2})$ , and $τ_{s} = \frac{1}{σ_{s}^{2}}$ . Meanwhile, temporal dynamics are modeled using a first order autoregressive process (AR(1)) with the following temporal relationship:

\begin{aligned} s_{i 1} = z_{i 1}, w h e r e s_{i t} = ϕ s_{i, t - 1} + z_{i t} w h e r e (t \geq 2) \end{aligned}

(10)

Where $ϕ \sim U (0.05, 0.95)$ .

The priors for the regression coefficients are presented in equation (11), where the coefficients are assumed to follow a normal distribution with mean 0 and standard deviation 5. Gelman³⁸ emphasize that a Normal (0, 5) prior is sufficiently broad to accommodate potential variations in effects, yet remains narrow enough to maintain Markov Chain Monte Carlo (MCMC) convergence stability and avoid the common issues associated with overly flat non-informative priors.

\begin{aligned} β_{k} & \sim N (0, 5^{2}), k = 1, \dots, p . \end{aligned}

(11)

\begin{aligned} σ_{δ} & \sim H a l f - N o r m a l (0, {2.5}^{2}), τ_{δ} = \frac{1}{σ_{δ}^{2}} \end{aligned}

(12)

The year effect is defined by setting $δ_{1}$ = 0, whereas for the other years:

\begin{aligned} δ_{t} \sim N (0, σ_{δ}^{2}), t = 2, \dots, T . \end{aligned}

(13)

The primary posterior quantity of interest is $θ_{i t}$ , which is defined on the log scale. The small area mean per capita expenditure on the original scale is obtained by exponentiating the posterior estimates of $θ_{i t}$ .

\begin{aligned} {\hat{μ}}_{i t} = \exp ({\hat{θ}}_{i t}) \end{aligned}

(14)

Parameter estimation was performed using the MCMC simulation approach with the Gibbs Sampling algorithm using R-Studio with the tipsae³⁶ and nimble.³⁹ The simulation was run for 10,000 iterations, with the first 5000 discarded as a burn-in period to ensure chain stability, and every 5^th iteration retained (thinning = 5) to reduce autocorrelation between samples.

3.4 Model evaluation

Model performance was evaluated using four primary indicators: the Coefficient of Variation (CV), Root Mean Squared Error (RMSE), Coverage Probability (CP), and a goodness of fit test based on the Wald statistic. Since the true population parameters are unobserved, model evaluation relies on the design-based direct estimates as a reference benchmark. Following Rao and Molina,¹¹ the CV, RMSE, and CP are calculated by treating the direct estimates as proxies for the true values.

The precision of the estimation results using the CV is calculated using the following formula:

\begin{aligned} C V = \frac{\sqrt{V ({\hat{μ}}_{i t})}}{{\hat{μ}}_{i t}} \times 100 % \end{aligned}

(15)

Where ${\hat{μ}}_{i t}$ represents the estimated average per capita expenditure and $V ({\hat{μ}}_{i t})$ denotes its variance. A smaller CV value indicates higher precision in the estimation.¹¹ Meanwhile, the estimation bias is measured using the RMSE according to the following formula:

\begin{aligned} R M S E = \sqrt{\frac{1}{m T} \sum_{i}^{m} \sum_{t}^{T} {({\hat{μ}}_{i t} - {\hat{Y}}_{i t}^{d i r})}^{2}} \end{aligned}

(16)

Where ${\hat{Y}}_{i t}^{d i r}$ denotes the design based direct estimate used as a proxy benchmark and ${\hat{μ}}_{i t}$ denotes the SAE estimate. The RMSE is used to assess how closely the estimated results approximate the actual values; a smaller RMSE indicates superior model accuracy.⁴⁰ Furthermore, model reliability is measured using the Coverage Probability (CP) at a 95% confidence level, defined as the proportion of confidence intervals that contain the actual parameter value. The CP formula is expressed in equation (17) as follows:

\begin{aligned} C P = \frac{1}{m T} \sum_{i = 1}^{m} \sum_{t = 1}^{T} I ({\hat{Y}}_{i t}^{d i r} \in C I_{i t}^{(1 - α)}) \end{aligned}

(17)

Where $C I_{i t}^{(1 - α)}$ denotes the $(1 - α)$ confidence interval. A desirable CP value should approximate the nominal level (e.g., 95%). An excessively low CP indicates under coverage, whereas an excessively high CP suggests that the intervals are overly broad.¹¹ Subsequently, an evaluation of model goodness of fit is conducted using the Wald statistic, calculated via the following formula:

\begin{aligned} W = \frac{1}{m T} \sum_{i = 1}^{m} \sum_{t = 1}^{T} \frac{{({\hat{μ}}_{i t} - {\hat{Y}}_{i t}^{d i r})}^{2}}{V ({\hat{μ}}_{i t}) - V ({\hat{Y}}_{i t}^{d i r})} \end{aligned}

(18)

Where ${\hat{μ}}_{i t}$ represents the estimated average per capita expenditure and $V ({\hat{μ}}_{i t})$ a is its variance, while ${\hat{Y}}_{i t}^{d i r}$ represents the direct estimate and $V ({\hat{Y}}_{i t}^{d i r})$ denotes the variance of the direct estimate. A smaller Wald value indicates a better model fit to the data.⁴¹

To strengthen the inferential assessment of model performance, we complement the descriptive metrics with a distribution-free comparison of predictive accuracy. Specifically, we apply a two sample Kolmogorov-Smirnov test to compare the distribution of absolute prediction errors between the best model and alternative models, without relying on parametric assumptions.⁴²

Table 2.

Pearson correlation analysis results between auxiliary variables and per capita expenditure.

Variable	Pearson correlation	p-value	Significance	Strength	Direction	Decision
NTL	0.503	0.000	Yes	Moderate	Positive	Yes
NDVI	−0.495	0.000	Yes	Moderate	Negative	Yes
NDBI	0.530	0.000	Yes	Moderate	Positive	Yes
MNDWI	0.351	0.000	Yes	Weak	Positive	No
LST	0.516	0.000	Yes	Moderate	Positive	Yes
CO	0.158	0.036	Yes	Very Weak	Positive	No
NO₂	0.424	0.000	Yes	Moderate	Positive	Yes
Distance to POI	−0.362	0.000	Yes	Weak	Negative	No
POI Density	0.492	0.000	Yes	Moderate	Positive	Yes

3.5 Benchmarking

Benchmarking was conducted using the raking method (iterative proportional fitting). This process is performed to meet official statistical standards, requiring that the final SAE estimates be consistent with aggregate figures published by BPS. The equation is as follows:

\begin{aligned} {\hat{Y}}_{i t}^{B E N C H} = {\hat{μ}}_{i t} + ({\hat{Y}}_{i t}^{d i r} - \sum_{i = 1}^{m} \sum_{t = 1}^{T} W_{i t} {\hat{μ}}_{i t} \end{aligned}

(19)

Where $W_{i t}$ represents the adjustment factor based on population size. The weight $W_{i t}$ is formulated as follows:

\begin{aligned} W_{i t} = \frac{N_{i t}}{\sum_{j ϵ g} N_{j t}} \end{aligned}

(20)

where

N_{i t}

represents the population share of subdistrict i within its corresponding regency/municipality in year t. The weights ensure that the population weighted aggregate of the benchmarked subdistrict estimates matches the official regency level estimate.

4 Results

The auxiliary variables used in the model are those exhibiting moderate to strong correlations (|r|>0.4) as correlation coefficients above 0.4 indicate a moderate relationship deemed sufficiently relevant for further analysis.³⁵ Table 2 shows that six variables NTL, NDVI, NDBI, LST, NO2, and POI Density were selected as auxiliary variables. Furthermore, variance smoothing was conducted. As illustrated in Figure 1, the difference is evident that prior to smoothing, the variance dispersion was extremely wide with extreme values; however, after smoothing, the variances were corrected to become more stable with a narrower dispersion.⁴³

This study compares four models: the HB model without spatio-temporal effects, the HB model with spatial effects, the HB model with temporal effects, and the HB model with spatio-temporal effects. The posterior mean results are presented in Table 3. Furthermore, Table 3 reveals that the relationship between the average posterior fixed-effect values and the average per capita expenditure varies across models. Notably, the model with spatio-temporal effects exhibits a pattern consistent with the prior correlation analysis. However, this observation alone does not demonstrate that the spatio-temporal model is superior to the others. Therefore, several diagnostic tests are necessary to assess model reliability and select the best model for benchmarking.

Figure 1.

Comparison of variance before smoothing (raw) and after smoothing (smoothed).

Table 3.

Posterior fixed effect estimation results of SAE model.

Model	HB non Spatio-Temporal	HB Spatial	HB temporal	HB spatio-Temporal
Variable	Mean	Mean	Mean	Mean
Intercept	13.254	14.232	13.981	12.793
NTL	−0.003	0.002	0.006	0.001
NDVI	−0.678	0.220	0.633	−0.199
NDBI	0.036	0.938	1.907	0.441
LST	0.017	−0.003	−0.000	0.020
NO₂	−0.049	−0.114	−0.080	0.074
POI Density	0.001	0.002	−0.000	0.000
rho	-	0.334	-	0.228
AR(1)			0.590	0.836

First, a check for model consistency is required, as the HB method employed applies computational sampling via MCMC. Figure 2 shows that each model is relatively consistent, with residuals tending to decrease as the sample size increases. In comparing models, the spatio-temporal model performs best, with residuals most tightly clustered around zero for both small and large samples. Next, Table 4 presents several indicators to assess and compare model performance which is complemented by an examination of the violin plot of the CV distribution in Figure 3. Furthermore, the results of the Kolmogorov–Smirnov tests (Table 5) provide inferential support for the observed performance differences. The p-values for all model comparisons are significantly lower than 0.05, indicating that the reduction in estimation error is statistically significant and consistent across the entire distribution. Taking all these indicators into account, it can be concluded that the HB Spatio-Temporal model is the strongest candidate. Overall, the proposed spatio-temporal framework demonstrates improved precision, reduced estimation bias, and a more coherent characterisation of uncertainty, supporting its suitability for small-domain estimation within official statistics.

Figure 2.

Scatterplot of residual model and number of samples.

Figure 3.

CV distribution using violinplot, boxplot and scatterplot.

Table 4.

Model evaluation.

Model	CV	RMSE	CP	Wald
Direct Estimate	13.21	-	-	-
HB Non Spatio-Temporal	5.84	290,118.0	64.20	4.72e + 15
HB Spatial	5.69	1,137,608.7	17.61	9.44e + 16
HB Temporal	2.18	723,702.9	12.50	2.52e + 17
HB Spatio-Temporal	9.31	224,845.7	97.16	1.83e + 15

Table 5.

Kolmogorov–Smirnov tests comparing distributions of absolute prediction errors.

Model comparison	K-S statistic	p-value
HB Spatio-Temporal vs HB Non Spatio-Temporal	0.4488	0.000
HB Spatio-Temporal vs HB Temporal	0.6477	0.000
HB Spatio-Temporal vs HB Spatial	0.7556	0.000

The HB Spatio-Temporal model was selected for further analysis. Consequently, benchmarking was conducted to ensure that the subdistrict level estimation results are consistent with the official aggregate results at the regency/municipality level. Table 6 presents the point estimates of the subdistrict average per capita expenditure before and after the benchmarking process. In general, it can be observed that the point estimates before and after benchmarking do not differ significantly. Specifically, the boxplot (left) shows that the distribution of per capita expenditure estimates across subdistricts remains relatively unchanged after benchmarking (Figure 4). The median values, quartile dispersion, and the presence of outliers are nearly identical between the pre and post benchmarking conditions. Furthermore, the scatterplot (right) also indicates a near perfect fit (R² = 1) between the two estimates.

Table 6.

Descriptive statistics of point estimates before and after benchmarking.

Benchmarking	Count	Min	Max	Mean	Median	SD
Before	176	700803.8	5254389.7	1515451.9	1315254.2	722877.1
After	176	686911.3	5250432.9	1518622.2	1317196.1	726253.7
Difference	176	−14956.6	56051.2	3170.4	−2368.3	21383.0

The consistency of the estimation results is further assessed through spatial visualisation. Figure 5 presents the spatial distribution of per capita expenditure across 44 subdistricts in Medan Municipality and Langkat Regency for the year 2024. The side-by-side comparison indicates that the spatial patterns remain highly consistent before and after the benchmarking process. Areas with higher expenditure levels continue to be concentrated in Medan, while predominantly subdistricts in Langkat remain associated with lower expenditure levels. These findings suggest that the benchmarking procedure introduces only limited local adjustments while preserving the underlying spatial relationships captured by the spatio-temporal hierarchical Bayesian small area estimation model.

Figure 4.

Comparison of estimates before and after benchmarking.

Figure 5.

Comparison of spatial distribution of subdistrict per capita expenditure in Medan and Langkat (2024): before benchmarking (left) and after benchmarking (right).

5 Discussion

This study demonstrates that applying Spatio-Temporal Hierarchical Bayesian Small Area Estimation (ST HB SAE) effectively improves the accuracy of average per capita expenditure estimates at the subdistrict level. Compared to direct estimation from Susenas, the spatio-temporal model provides significantly more precise and stable results, particularly in subdistricts with small sample sizes. These findings confirm that modeling based on “borrowing strength” through spatial (CAR) and temporal (AR1) components can overcome the primary limitations of official surveys in small domains.

The integration of geospatial big data has also proven to provide a significant contribution to strengthening the SAE model. Variables such as NTL, NDBI, NDVI, LST, and NO₂ exhibit consistent correlations with household expenditure, thereby enriching the model's information. These indicators capture variations in physical, environmental, and economic activity not covered in official survey data or conventional administrative records. This result aligns with recent studies emphasizing the potential of remote sensing indicators for measuring socioeconomic dynamics and demonstrates that integrating geospatial big data offers tangible added value in the Indonesian context, which features contrasting urban rural characteristics.

In the context of model implementation, ST HB SAE exhibits superior performance compared to nonspatial models and spatial models without a temporal component. The improvement in precision is reflected in the significant reductions in CV and RMSE, as well as in the high coverage probability. Furthermore, the consistency of results against official aggregates following the benchmarking process (adjustment < 1 percent) indicates that this model aligns with official statistical standards and can be integrated into existing statistical production frameworks.

From a regional substantive perspective, spatial patterns reveal a higher concentration of per capita expenditure in Medan's urban centers, whereas agrarian areas in Langkat consistently remain in the lower quantiles. Temporal trends from 2021 to 2024 show a stable increase in household expenditure, indicating the post-pandemic economic recovery. The spatio-temporal model successfully captures these patterns more smoothly and realistically compared to direct estimation, which tends to be fluctuating.

Overall, the findings of this study support the initial research objectives that the ST HB SAE model produces average per capita expenditure estimates that are more precise, temporally stable, and sensitive to spatial differences. From a policy perspective, this approach offers a stronger foundation for local governments and BPS to monitor welfare dynamics at the subdistrict level, strengthen the basis for targeting social programs, and drive the modernization of official statistics through the utilization of geospatial big data and Bayesian methodologies. Nevertheless, this study also offers opportunities for further development, such as integrating higher spatial resolution, using dynamic administrative data, or expanding the model to other regions to assess the method's generalizability.

6 Conclusions

This study demonstrates that applying ST HB SAE effectively improves the accuracy of average per capita expenditure estimates at the subdistrict level in Medan Municipality and Langkat Regency. Compared to Susenas direct estimates, the spatio-temporal model produces estimates that are more precise, temporally stable, and capable of covering non-sampled subdistricts. The geospatial auxiliary variables used in this study NTL, NDVI, NDBI, LST, NO₂, and POI indicators show statistically significant relationships with household expenditure, as indicated by correlation tests. Although the correlation strengths vary across variables, all indicators contribute positively to the model's predictive capability. These results underscore the potential of geospatial data as a relevant source of supplementary information for small area estimation, particularly in regions lacking annually updated administrative data.

The application of ST HB SAE also yields superior performance compared to nonspatial models and spatial models without a temporal component, as reflected in reduced CV and RMSE and increased coverage probability. The consistency with official aggregates following benchmarking demonstrates that this model meets the requirements of official statistics and can serve as a complement to welfare estimation at the subdistrict level.

Overall, these findings support the use of the Susenas panel, combined with geospatial big data, to produce more accurate small area estimates. This model has the potential to assist local governments in identifying regional inequalities more precisely and strengthening evidence based development planning. Going forward, future research could explore integrating higher resolution geospatial data, incorporating dynamic administrative variables, or applying the model to other regions to test the generalizability of this approach.

Footnotes

Acknowledgments

This manuscript is a portion of the first author's doctoral research. The authors wish to express their gratitude to The Graduate School, Universitas Gadjah Mada, for granting access to literature and research resources, to BPS for providing the data, and to the Indonesia Endowment Fund for Education or Lembaga Pengelola Dana Pendidikan (LPDP) under the Ministry of Finance of the Republic of Indonesia for sponsoring the author's doctoral degree.

ORCID iDs

Hasrat Ifolala Zebua

Sigit Heru Murti Budi Santosa

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Zebua

. Pemodelan Kemiskinan Di Sumatera Utara Menggunakan Regresi Nonparametrik Kernel Dan Splines (Kernel and splines nonparametric regression for poverty model in North Sumatera). Proc Semnas Offstat 2021; 2021: 899–907.

World Bank. Poverty and shared prosperity 2022. Washington, DC: Correcting Course, 2022.

UNECE. Poverty measurement. Geneva, Switzerland: Guide to Data Disaggregation, 2020.

BPS. Profil Kemiskinan Di Indonesia Maret 2024. Jakarta, Indonesia: Badan Pusat Statistik RI, 2024.

Suparman

Muzakir

. Regional inequality, human capital, unemployment, and economic growth in Indonesia: Panel regression approach. Cogent Economics & Finance 2023; 11: 2251803.

Kanbur

Venables

. Spatial inequality and development. Newyork: Oxford University Press, 2005.

Sachs

McArthur

Schmidt-Traub

, et al. Ending Africa’s poverty trap. Brookings pap. Econ Act 2004; 2004: 117–240.

BPS. Data Dan Informasi Kemiskinan Kabupaten/Kota Tahun 2023. vol. 15. Jakarta: Badan Pusat Statistik, 2023.

BPS. Penghitungan Dan Analisis Kemiskinan Makro Di Indonesia Tahun 2024. (2024).

10.

Pratesi

. Analysis of poverty data by small area estimation. United Kingdom: John Wiley & Sons Ltd, 2016.

11.

Rao

JNK

Molina

. Small area estimation: second edition. 2015. doi:10.1002/9781118735855

12.

Masaki

Newhouse

Silwal

, et al. Small area estimation of non-monetary poverty with geospatial data. 2020. doi:10.1596/1813-9450-9383

13.

Aziz

Ubaidillah

. Big data for small area estimation: happiness index with twitter data. Proc Int Conf Data Sci Off Stat 2021; 2021: 963–972.

14.

NASEM. Toward a 21st century national data infrastructure: enhancing survey programs by using multiple data sources. Washington, D.C: National Academies Press, 2023, doi:10.17226/26804.

15.

Kaban

Nasution

Caraka

, et al. Implementing night light data as auxiliary variable of small area estimation. Commun. Stat. Theory Methods 2024; 53: 310–327.

16.

Putri

Wijayanto

Sakti

. Developing relative spatial poverty index using integrated remote sensing and geospatial big data approach: a case study of east Java, Indonesia. ISPRS Int J Geoinf 2022; 11: 275.

17.

Putri

Wijayanto

Pramana

. Multi-source satellite imagery and point of interest data for poverty mapping in East Java, Indonesia: Machine learning and deep learning approaches. Remote Sens Appl 2023; 29: 100889.

18.

Utami

Wijayanto

Pramana

. Spatially granular poverty index (SGPI) for urban poverty mapping in Jakarta metropolitan area (JMA): a remote sensing satellite imageries and geospatial big data approach. Earth Sci Inf 2023; 16: 3531–3544. Preprint at https://doi.org/10.1007/s12145-023-01084-7.

19.

Permatasari

Laksono

Ubaidillah

. Small area estimation of poverty using remote sensing data (Case Study: Expenditure Per Capita Estimation of Very Poor Households in West Java, Indonesia). In: 64th ISI world statistics congress, 2023.

20.

Tian

Yang

, et al. Highway networks and regional poverty: evidence from Chinese counties. Structural Change and Economic Dynamics 2024; 69: 224–231.

21.

Hoogstra

Velichety

Zhang

. Developing a contextual model of poverty prediction using data science and analytics – the case of Shelby County. Decis Support Syst 2024; 177: 114080.

22.

Tang

Zhao

Zhang

, et al. Poverty estimation using a ConvLSTM-based model with multisource remote sensing data: a case study in Nigeria. IEEE J Sel Top Appl Earth Obs Remote Sens 2024; 17: 3516–3529.

23.

Arshad

Zulfiqar

Zaib

, et al. Mapping socioeconomic conditions using satellite imagery: a computer vision approach for developing countries. Journal of Economy and Technology 2023; 1: 144–163.

24.

Han

Zhang

Zhao

, et al. High-resolution remote sensing data can predict household poverty in pastoral areas, Inner Mongolia. China Geography and sustainability 2021; 2: 254–263.

25.

Subash

Kumar

Aditya

. Satellite data and machine learning tools for predicting poverty in rural India. Agricultural Economics Research Review 2018; 31: 231–240.

26.

Jean

, et al. Combining satellite imagery and machine learning to predict poverty. http://science.sciencemag.org/ (2016).

27.

Duque

Patino

Ruiz

, et al. Measuring intra-urban poverty using land cover and texture metrics derived from remote sensing data. Landsc Urban Plan 2015; 135: 11–21.

28.

Xie

Jean

Burke

, et al. Transfer learning from deep features for remote sensing and poverty mapping. Proceedings of the AAAI Conference on Artificial Intelligence 2016; 30: 3929–3935.

29.

Lee

Braithwaite

. High-resolution poverty maps in Sub-Saharan Africa. World Dev 2022; 159: 106028.

30.

Tan

. Estimation of poverty based on remote sensing image and convolutional neural network. Advances in Remote Sensing 2019; 08: 89–98.

31.

Marchetti

, et al. Small area model-based estimators using big data sources. J Off Stat 2015; 31: 263–281.

32.

Gujarati

. Basic econometrics, 5th ed. New York, NY: McGraw-Hill, 2009.

33.

Sugiyono. Statistika Untuk Penelitian. Bandung: Alfabeta, 2007.

34.

Mukaka

. Statistics corner: a guide to appropriate use of correlation coefficient in medical research. Malawi Med J 2012; 24: 69–71.

35.

Evans

. Straightforward statistics for the behavioral sciences. Pacific Grove, CA: Thomson Brooks/Cole Publishing Co, 1996.

36.

De Nicolò

Gardini

. The R package Tipsae: Tools for mapping proportions and indicators on the unit interval. J Stat Softw 2024; 108. doi:10.18637/jss.v108.i01

37.

White

McConville

Moisen

, et al. Hierarchical Bayesian small area estimation using weakly informative priors in ecologically homogeneous areas of the interior western forests. Frontiers in Forests and Global Change 2021; 4: 752911.

38.

Gelman

Jakulin

Pittau

, et al. A weakly informative default prior distribution for logistic and other regression models. Ann Appl Stat 2008; 2: 1360–1383.

39.

NIMBLE Development Team. NIMBLE: MCMC, Particle Filtering, and Programmable Hierarchical Modeling. Preprint at (2024).

40.

Chai

Draxler

. Root mean square error (RMSE) or mean absolute error (MAE)? – arguments against avoiding RMSE in the literature. Geosci Model Dev 2014; 7: 1247–1250.

41.

McCulloch, Charles

Searle, Shayle

Neuhaus

. Generalized, linear, and mixed models, second edition. New York: John Wiley & Sons, 2001.

42.

Hassani

Silva

. A kolmogorov-smirnov based test for comparing the predictive accuracy of two sets of forecasts. Econometrics 2015; 3: 590–609.

43.

De Nicolò

Gardini

. The R package Tipsae: tools for mapping proportions and indicators on the unit interval [Computer software manual] 2024.