A Random Forest–Augmented Fay–Herriot Model for Nonlinear Area-Level Small Area Estimation

2.1 Conceptual basis of the proposed Random Forest–augmented Fay–Herriot model

Machine learning has recently been integrated into SAE, most notably through the MERF of Krennmair and Schmid [12], which extends the mixed-model framework of Hajjem et al. [23]. MERF modifies the EM algorithm of Wu and Zhang [24] by replacing the fixed-effect term $X_{i}^{T}\beta $ with OOB Random Forest predictions. While MERF works effectively for unit-level SAE, it cannot be directly transferred to the area-level Fay–Herriot model due to two reasons: (1) Fay–Herriot responses are area means, not unit-level outcomes; and (2) Random Forest minimizes prediction error [25], whereas the EM algorithm maximizes likelihood; thus, their objectives are not naturally aligned.

To retain the likelihood principles of Fay–Herriot while enabling nonlinear prediction, the proposed approach applied Random Forest after the FH–EM algorithm converges. This ensures that variance estimation follows the standard Fay–Herriot likelihood and the fixed-effect predictor is either replaced or corrected using Random Forest, depending on the variant considered.

The original Fay–Herriot model [26] is given in Eq. (1):

${{\hat{\theta }}_{i}}=X_{i}^{T}\beta +{{z}_{i}}{{v}_{i}}+{{e}_{i}}$,$~i=1,\cdots ,m$       (1)

where, $i$ represents the area; ${{X}_{i}}$ is a $p\times 1$ vector of area-level covariates; $\beta $ is the corresponding vector of regression parameters; $v_i \sim N\left(0, \sigma_v^2\right)$ is a random component which captures the unobserved heterogeneity among areas; and $e_i \sim N\left(0, \sigma_{e_i}^2\right)$ is sampling error, where $\sigma _{{{e}_{i}}}^{2}$ is known from the sampling design; and ${{z}_{i}}=1$ in the classical Fay–Herriot specification. The covariance structure is $\Omega =ZG{{Z}^{T}}+S$, with $G={{I}_{m}}\sigma _{v}^{2}$ and $S={{I}_{m}}\sigma _{e}^{2}$, where ${{I}_{m}}$ denotes the $m\times m$ identity matrix.

The key idea of the FH–RF approach is to retain the random-effects structure as expressed in Eq. (2):

${{\hat{\theta }}_{i}}=(X_{i}^{T}\beta )+\left( {{z}_{i}}{{v}_{i}}+{{e}_{i}} \right)$         (2)

While substituting or correcting the linear fixed-effect component with a Random Forest predictor to capture nonlinear relationships that may be present in socioeconomic and area-level data [27].

2.2 Expectation–maximization algorithm for the Fay–Herriot model

Valdez et al. [22] proposed a Fay–Herriot-specific EM algorithm to ensure stable estimation of $\sigma _{v}^{2}$ and avoid zero/negative values. The algorithm iterations are shown in Eqs. (3)–(5):

• E-step (Eq. (3)):

${{\hat{v}}^{\left( r+1 \right)}}={{\left( \frac{1}{\hat{\sigma }_{v}^{2,\left( r \right)}}{{I}_{m}}+{{S}^{-1}} \right)}^{-1}}{{S}^{-1}}\left( \hat{\theta }-{{X}^{T}}{{\beta }^{\left( r \right)}} \right)$         (3)

• M-step (Eq. (4) and Eq. (5)):

${{\hat{v}}^{\left( r+1 \right)}}={{\left( \frac{1}{\hat{\sigma }_{v}^{2,\left( r \right)}}{{I}_{m}}+{{S}^{-1}} \right)}^{-1}}{{S}^{-1}}\left( \hat{\theta }-{{X}^{T}}{{\beta }^{\left( r \right)}} \right)$       (4)

$\hat{\sigma }_{v}^{2,\left( r+1 \right)}=\frac{1}{m}\left( {{{\hat{v}}}^{\left( r+1 \right)T}}{{{\hat{v}}}^{\left( r+1 \right)}}+tr{{\left( \frac{1}{\hat{\sigma }_{v}^{2,\left( r \right)}}{{I}_{m}}+{{S}^{-1}} \right)}^{-1}} \right)$           (5)

The updating process is repeated until the change in the parameter estimates becomes very small, for instance, below 10⁻⁵. The EM output is then used as the foundation for integrating Random Forest to produce the two model variants below.

2.3 Proposed model 1: Fay–Herriot with Random Forest correction

FHRF-C retains the linear Fay–Herriot predictor but augments it with a nonlinear Random Forest correction, as shown in Eq. (6):

${{\hat{\theta }}_{i}}={{f}^{*}}\left( {{X}_{i}} \right)+{{v}_{i}}+{{e}_{i}}$         (6)

where, the corrected fixed-effect component is defined in Eq. (7):

${{f}^{*}}\left( {{X}_{i}} \right)=X_{i}^{T}\hat{\beta }+\hat{f}\left( {{X}_{i}} \right)$        (7)

Steps:

1. Perform FH–EM to obtain $\hat{\beta }$, ${{\hat{v}}_{i}}$, and $\hat{\sigma }_{v}^{2}$.
2. Construct residuals using Eq. (8):

${{\hat{\varepsilon }}_{i}}={{\hat{\theta }}_{i}}-X_{i}^{T}\hat{\beta }-{{\hat{v}}_{i}}$,  $i=1,\cdots ,m$           （8）

3. Fit Random Forest to$~\left( {{{\hat{\varepsilon }}}_{i}},~{{X}_{\left( i \right)}} \right)~$using bootstrap sampling (SRSWR) to construct the Random Forest model. Predictions $\hat{f}\left( {{X}_{\left( i \right)}} \right)$, as shown in Eq. (9), rely solely on OOB estimates to avoid overfitting:

$\hat{f}\left( {{X}_{i}} \right)={{\left[ \hat{f}\left( {{X}_{1}} \right),~\cdots ,\hat{f}\left( {{X}_{m}} \right) \right]}^{T}}$       (9)

4. The corrected fixed-effect ${{f}^{*}}\left( {{X}_{i}} \right)$ is then computed using Eq. (7).
5. SAE is carried out using the Empirical Best Linear Unbiased Predictor (EBLUP) [1] in Eq. (10):

$\begin{aligned} \hat{\theta}_i^{F H R F-C} & =\left(T\left(\hat{\theta}_i\right)\right)\left(\sigma_v^2\right) =f^*\left(\boldsymbol{X}_{\boldsymbol{i}}\right)+\frac{\sigma_v^2}{\left(\sigma_v^2+\sigma_{e_i}^2\right)}\left(\hat{\theta}_i-f^*\left(\boldsymbol{X}_{\boldsymbol{i}}\right)\right)\end{aligned}$       (10)

This model is suitable when the Fay–Herriot linear predictor is moderately misspecified but still informative.

2.4 Proposed model 2: Fay–Herriot with Random Forest replacement

FHRF-R replaces the entire Fay–Herriot linear predictor with a nonparametric Random Forest function, yielding Eq. (11):

${{\hat{\theta }}_{i}}=\hat{f}\left( {{X}_{i}} \right)+{{v}_{i}}+{{e}_{i}}$       (11)

Steps:

1. Perform FH–EM to obtain $\hat{\beta }$, ${{\hat{v}}_{i}}$, and $\hat{\sigma }_{v}^{2}$.
2. Construct an adjusted response using Eq. (12):

$\hat{\theta }_{i}^{*}={{\hat{\theta }}_{i}}-{{\hat{v}}_{i}}$, $i=1,\cdots ,m$       (12)

3. Fit Random Forest to $\left( \hat{\theta }_{i}^{*},~{{X}_{\left( i \right)}} \right)$ using OOB predictions.
4. Obtain Random Forest fixed-effect component $\hat{f}\left( {{X}_{i}} \right)$ as in Eq. (13):

$\hat{f}\left( {{X}_{i}} \right)={{\left[ \hat{f}\left( {{X}_{1}} \right),~\cdots ,\hat{f}\left( {{X}_{m}} \right) \right]}^{T}}$          (13)

5. Computation of small area estimates using the EBLUP as shown in Eq. (14):

$\begin{aligned} \hat{\theta}_i^{F H R F-R} & =\left(T\left(\hat{\theta}_i\right)\right)\left(\sigma_v^2\right)=\hat{f}\left(\boldsymbol{X}_{\boldsymbol{i}}\right)+\frac{\sigma_v^2}{\left(\sigma_v^2+\sigma_{e_i}^2\right)}\left(\hat{\theta}_i-\hat{f}\left(\boldsymbol{X}_{\boldsymbol{i}}\right)\right)\end{aligned}$       (14)

This variant is appropriate when the Fay–Herriot linearity assumption is strongly violated.

2.5 Practical selection between Fay–Herriot with Random Forest correction and Fay–Herriot with Random Forest replacement

Although FHRF-C and FHRF-R differ conceptually in how Random Forest is incorporated into the Fay–Herriot framework, practical applications require guidance on model selection. In general, FHRF-C is preferable when the relationship between area-level covariates and the target parameter is predominantly linear with localized nonlinear residual patterns, allowing the classical Fay–Herriot predictor to remain effective. The Random Forest component captures departures from linearity. In contrast, FHRF-R is more suitable when the mean structure itself is strongly nonlinear, so that replacing the linear predictor with a full Random Forest function provides a better fit. Diagnostic evidence such as significant RESET or Generalized Additive Model (GAM) test results, pronounced systematic residual patterns, or substantial reductions in Random Forest OOB error relative to the linear Fay–Herriot model may indicate preference for FHRF-R, whereas milder nonlinear residual structures favor FHRF-C. From a computational perspective, FHRF-C is typically more efficient because it retains the closed-form Fay–Herriot predictor, while FHRF-R requires repeated Random Forest estimation. A summary of practical diagnostic indicators and trade-offs is provided in Table 1.

Table 1. Summarizes practical diagnostic indicators and trade-offs for selecting between Fay–Herriot with Random Forest correction (FHRF-C) and Fay–Herriot with Random Forest replacement (FHRF-R)

2.6 Bootstrap mean squared error estimation for Random Forest–Augmented Fay–Herriot models

Analytical MSE formulas are not available once Random Forest components replace or correct the linear predictor. Building on the methodology proposed by Rao and Molina [1], the bootstrap method is employed to evaluate the MSE of the FHRF models. The main steps of this procedure are summarized below:

1. Fit FHRF to obtain the fixed component $\hat{f}\left( {{X}_{i}} \right)$ and $\hat{\sigma }_{v}^{2}$.
2. For b = 1, …, B:

1. Generate pseudo-true area values using (15):

$\theta_{i *}^{(b)} \sim N\left(\hat{f}\left(\boldsymbol{X}_{\boldsymbol{i}}\right), \hat{\sigma}_v^2\right)$      (15)

2. Generate synthetic direct estimates using (16):

$\hat{\theta}_{i *}^{(b)} \sim N\left(\theta_{i *}^{(b)}, \sigma_{e_i}^2\right)$        (16)

3. Refit FHRF and obtain $\hat{\theta }_{i*}^{{{H}^{\left( b \right)}}}$.

3. Compute bootstrap MSE using Eq. (17):

$\widehat{mse}\left( {{\theta }_{i}} \right)=ms{{e}_{B}}\left( {{{\hat{\theta }}}_{i}} \right)=\frac{1}{B}\underset{b=1}{\overset{B}{\mathop \sum }}\,{{\left( \hat{\theta }_{i*}^{{{H}^{\left( b \right)}}}-\theta _{i*}^{\left( b \right)} \right)}^{2}}$        (17)

4. Compute the relative standard error (RSE) using Eq. (18):

$\widehat{RSE}\left( {{\theta }_{i}} \right)=\frac{\sqrt{\widehat{mse}\left( {{\theta }_{i}} \right)}}{\hat{\theta }_{i*}^{{{H}^{\left( b \right)}}}}\times 100%$         (18)

Bootstrap MSE captures uncertainty from both the sampling design and the Random Forest-enhanced model structure.

4.1 Simulation results

Algorithmic convergence is an essential consideration in EM-based SAE models [28]. To ensure stable initialization, the proposed method adopts starting values from an EBLUP-based Fay–Herriot formulation [1], which is computationally efficient and structurally consistent with the FH–EM framework. The simulation consisted of 100 replications, each covering population generation, sampling, direct estimation, and model fitting. Performance is evaluated using RRMSE (precision), RB (bias), and replication-based stability measures MRAE and MdRAE.

Table 3 presents the RRMSE values across all simulation settings (L1–L5, NL1–NL5). Both FHRF variants consistently dominate direct estimation, achieving uniformly lower mean and median RRMSE. Because the RRMSE differences violated normality (Shapiro–Wilk), the Wilcoxon Signed-Rank test [29] was employed. All comparisons yield extremely small p-values (7.19 × 10⁻³³ – 5.25 × 10⁻⁴³), confirming the statistical significance of the improvements. Beyond statistical significance, the magnitude of the gains is also substantial: the median RRMSE reduction from direct estimation to FHRF-R/FHRF-C ranged from 47.4% to 63.9% across scenarios. The largest reduction (63.9%) occurred in scenario L4, where the RRMSE decreased from 0.5357 under direct estimation to 0.2069 for the FHRF-C model, indicating a practically meaningful improvement in estimation accuracy.

Under nonlinear scenarios, FHRF-C and FHRF-R outperform FH–EM, as indicated by consistently lower mean and median RRMSE across all evaluated settings (Table 3). Specifically, both models achieved average RRMSE reductions of 11%–21% and median reductions of 14%–27% relative to the FH–EM estimator. Pairwise Wilcoxon signed-rank tests comparing FHRF-C and FHRF-R with FH–EM confirm that these improvements are statistically significant, with p-values ranging from 1.34 × 10⁻³ to 8.78 × 10⁻¹³. Beyond statistical significance, the magnitude of the improvements is also meaningful: the median RRMSE reduction from FH–EM to the FHRF models ranged from 14.3% to 27.2% across nonlinear scenarios. The largest reduction (27.2%) occurred in scenario NL5, where the RRMSE decreased from 0.2192 under FH–EM to 0.1595 for the FHRF-R model, highlighting the practical advantage of Random Forest in capturing nonlinear relationships without requiring explicit functional form specification [14, 15, 22].

Table 4 and Table 5 further examine robustness when the number of areas (m) and the number of covariates (p) are varied while preserving the same nonlinear data-generating mechanism (scenario NL1). Increasing the number of areas substantially strengthens the relative efficiency of the FHRF estimators: the average RRMSE reduction relative to FH–EM rises from about 16% at m = 100 to approximately 30–40% at m = 500, indicating that the machine-learning component benefits from richer cross-area information and more stable learning of nonlinear relationships.

Table 3. Relative root mean squared error (RRMSE) results across 100 simulated replications under linear and nonlinear scenarios

Note: SD = standard deviation. The smallest are marked in bold; FH–EM: Fay–Herriot model estimated via the expectation–maximization algorithm; FHRF-C: Fay–Herriot with Random Forest correction; FHRF-R: Fay–Herriot with Random Forest replacement.

Table 4. Relative root mean squared error (RRMSE) across 100 simulation replications under varying numbers of areas and covariates

Note: FH–EM: Fay–Herriot model estimated via the expectation–maximization algorithm; FHRF-C: Fay–Herriot with Random Forest correction; FHRF-R: Fay–Herriot with Random Forest replacement.

Table 5. Percentage reduction in relative root mean squared error (RRMSE) relative to FH–EM under varying dimensional settings

Note: FH–EM: Fay–Herriot model estimated via the expectation–maximization algorithm; FHRF-C: Fay–Herriot with Random Forest correction; FHRF-R: Fay–Herriot with Random Forest replacement.

Table 6. Summarizes of replication-level stability metrics for each nonlinear scenario, where the reported values represent the means across simulation replications

Note: MRAE: mean relative absolute error; MdRAE: median relative absolute error; FH–EM: Fay–Herriot model estimated via the expectation–maximization algorithm; FHRF-C: Fay–Herriot with Random Forest correction; FHRF-R: Fay–Herriot with Random Forest replacement.

tu_pian_28.png

Figure 1. Mean relative absolute error (MRAE) across 100 simulation replications (250 areas each) under the nonlinear scenario group

tu_pian_29.png

In linear scenarios, FH–EM remains the best-performing model, producing the lowest RRMSE values across L1–L5 (Table 3), indicating that FH–EM remained the best-performing estimator when the linear mixed-model assumptions were correctly specified. This result was theoretically expected because FH–EM is a maximum likelihood estimator under the Gaussian–linear mixed model, thereby enjoying optimal properties such as consistency and efficiency. This pattern was consistent with the results of Krennmair and Schmid [12], who found that unit-level SAE models incorporating Random Forest (SAE-MERF) outperform classical SAE only when the underlying relationships are nonlinear. When the response–covariate relationship is truly linear, the inclusion of nonparametric learners such as Random Forest does not provide additional benefits. It may have introduced approximation noise due to data-driven partitioning and bootstrap aggregation. Hence, an important limitation of the proposed FHRF approach was that it was not designed to improve performance under correctly specified linear settings, but rather to address cases where linear specifications are inadequate.

A direct comparison between FHRF-C and FHRF-R shows that both models produced similar RRMSE distributions and consistently outperformed the benchmarks in all nonlinear scenarios (NL1–NL5). Although Table 3 reports marginally lower mean or median RRMSE for FHRF-R in a few scenarios (NL1, NL2, NL4, NL5), the Wilcoxon Signed-Rank tests indicated that FHRF-C achieved significantly lower RRMSE in almost all cases, except NL2, where the two models are statistically indistinguishable. In terms of effect size, the differences between the two variants are generally small: the median RRMSE of FHRF-R is lower than that of FHRF-C by approximately 0.13%–0.94% in scenarios L3, NL1, NL2, and NL4, whereas FHRF-C shows larger advantages of about 1.13%–7.23% in scenarios L1, L2, L4, L5, NL3, and NL5. Replication-level error patterns (Figure 1) further show that both variants displayed comparable stability in terms of MRAE. Overall, these results suggested that while FHRF-C exhibits a modest statistical advantage, both FHRF variants remain equally dependable under nonlinear conditions.

While RRMSE captures both estimator accuracy and precision, examining RB remains essential for understanding systematic deviation. Interpreting RB alongside RRMSE reflects the classical bias–variance trade-off, in which an effective estimator must control both systematic bias and sampling variability. The RB patterns in Figure 2 indicate that direct estimates have a median close to zero, with a narrow interquartile range (IQR), because positive and negative sampling errors tend to offset each other under SRS [30]. However, the relatively large RRMSE of direct estimates reflects their high sampling variability, making direct estimation inefficient for small or heterogeneous domains despite their near-zero RB.

Under linear scenarios (L1–L5), FH–EM exhibits the most favorable RB, with medians near zero and smaller IQRs than both FHRF variants (Figures 2(a)–(e)), confirming its suitability when the actual relationship is linear. In contrast, under nonlinear scenarios (NL1–NL5), FH–EM exhibited systematic underestimation (medians below zero) and a wider IQR (Figures 2(f)–(j)), indicating difficulty in capturing nonlinear structure. Both Random Forest-based models addressed this limitation: FHRF-C and FHRF-R yielded RB distributions that are more tightly concentrated around zero, with IQRs reduced by approximately 41%–59% relative to FH–EM across all nonlinear scenarios, and exhibit nearly identical bias patterns.

Taken together, these results showed that the relative performance of FH–EM and FHRF models is fundamentally driven by the underlying functional form. FH–EM remained the preferred estimator under linear data-generating processes, whereas FHRF-C and FHRF-R provided clear advantages only when nonlinear relationships are present. Accordingly, the proposed FHRF framework should be viewed not as a universal replacement for the classical Fay–Herriot model, but as a complementary extension designed for nonlinear and heterogeneous settings frequently encountered in socioeconomic applications.

4.2 Analysis of empirical data

Building on the simulation evidence, the proposed methods were applied to an empirical case study of agricultural household expenditure in Indonesia. At the regency/city level, Statistics Indonesia commonly publishes average per capita expenditure estimates based on direct estimation [31]. While this approach is adequate for the general population, where RSE values typically fall within acceptable limits [1, 26], it becomes unreliable for subpopulations, such as agricultural households, where sample sizes are substantially smaller. Consistent with standard publication criteria (RSE ≤ 25%) [32], 15 regencies/cities exceeded this threshold, indicating the necessity of model-based SAE methods [28, 29, 33, 34].

Several diagnostic issues further confirmed the limitations of classical linear SAE approaches. First, clear nonlinear patterns were observed in the association between the outcome variable and the auxiliary covariates: the residuals failed the Shapiro–Wilk normality test, and both the GAM and Ramsey’s RESET tests rejected linear functional form assumptions (p < 0.05). Second, severe multicollinearity was observed among the auxiliary variables [35]. Third, the response variable was highly skewed (Z-skew = 35.548) [36]. Fourth, the unavailability of recent unit-level census covariates [37], due to Indonesia’s 10-year census cycle, restricted the feasibility of unit-level SAE. These issues collectively indicated that (i) direct estimation is imprecise for small domains, (ii) FH–EM is mis-specified for nonlinear structures, and (iii) area-level, machine-learning–enhanced models are operationally preferred.

Given these conditions, Random Forest offered several advantages: it captures nonlinear patterns [25], mitigates multicollinearity through randomized feature selection [38], and remains robust to skewness and non-normality [21, 33, 34]. Combined with the Fay–Herriot framework [26], these properties motivate the FHRF approach for this empirical context.

Both FHRF-C and FHRF-R yielded similar average RSE levels, at 6.56% (median 5.37%) and 6.35% (median 5.35%), respectively. The normality of paired differences was rejected by the Shapiro–Wilk test, necessitating the use of the Wilcoxon signed-rank test, which did not indicate a significant difference in median RSE (p = 0.144). In terms of effect size, the difference between the two models is small: FHRF-R achieved slightly lower RSE values, with a reduction of about 0.02 percentage points in the median (approximately 0.5% relative reduction) and 0.21 percentage points in the mean (approximately 3.2% relative reduction), indicating practically similar performance between the two approaches. Despite this, visual inspection showed that FHRF-R tended to produce more uniform gains, particularly among areas with initially high RSE. Most direct estimates fell below the 45° line when compared with FHRF-C (Figure 3(a)), indicating lower RSEs for FHRF-C, although a small number of points lay above the line. In contrast, almost all high-RSE direct estimates fell below the 45° line when compared with FHRF-R (Figure 3(b)), highlighting the superior performance of FHRF-R in these cases. Figure 3(c) shows that FH–EM produced higher RSEs than the direct estimator for many areas. The comparison in Figure 3(d) indicates that FHRF-R generally achieved lower RSEs than FHRF-C, particularly in areas with high sampling variability. Comparisons with FH–EM (Figures 3(e) and 3(f)) further show that both FHRF-C and FHRF-R generally outperformed FH–EM, with the advantage being more pronounced for FHRF-R.

tu_pian_30.png

The divergence between empirical and simulated findings may result from the stronger nonlinear structure in the real data and the concentration of problematic areas among high-RSE domains—features that were not explicitly emphasized in the simulations. Because the main empirical objective is to reduce the proportion of areas exceeding the 25% RSE threshold, subsequent analysis focuses on FHRF-R.

The SAE-FHRF-R model delivered clear gains in precision. The average RSE decreased from 7.46% (median 5.57%) under direct estimation to 6.35% (median 5.35%) with FHRF-R, corresponding to an approximate 15% reduction. The statistical significance of this reduction was confirmed by the Wilcoxon signed-rank test (p = 2.93 × 10⁻³⁸). In terms of effect size, this corresponds to an absolute reduction of 1.11 percentage points in the mean RSE (approximately 14.9% relative reduction) and 0.22 percentage points in the median (approximately 4.0% relative reduction), indicating a practically meaningful improvement in precision. The association between direct and model-based RSE is also strong (Pearson r = 0.852). This pattern is visible in Figure 3b, where most points fall below the 45° benchmark, reflecting systematic improvement. The most significant gains occurred in areas with high direct-estimate RSE (15%–60%), where the RSE values were consistently reduced to below the 25% publication threshold. In this context, the number of areas exceeding the 25% threshold decreases from 15 under direct estimation to 3 under FHRF-R. In contrast, areas with already low RSE exhibited minimal change. This selective improvement reflects the model’s adaptive behavior: As implied by Eq. (14), the fixed-effect component exerts greater influence when sampling-error variance (and thus RSE) is large, yielding marked reductions where they are most needed.

Comparisons with FH–EM further underscore this advantage. FH–EM produced a higher average RSE (6.80%) than FHRF-R (6.35%), and several areas exhibited RSE values even larger than their direct estimates (Figure 3c), indicating a lack of model fit. Across the full distribution, FHRF-R consistently outperformed FH–EM (Figure 3f), aligning with evidence that machine-learning methods are more effective at capturing complex and nonlinear relationships [25].

When areas are sorted by sample size, high RSE values from direct estimation—often above the 25% publication threshold—were concentrated in areas with very small agricultural household samples (Figure 4). In these domains, direct estimates became unstable, whereas FHRF-R consistently reduced RSEs to more acceptable levels. In such high-variability domains with small effective sample sizes, the reduction in RSE achieved by the FHRF-R model can reach up to 76%.

tu_pian_31.png

Figure 5 highlights that the most significant discrepancies between direct and FHRF-R estimates occur in areas with unusually high per capita expenditure relative to their neighbors. A clear example is Kota Manado (7171): the direct estimate was IDR 2,206,176 with an RSE of 34.07%, far higher than other districts in the province (IDR 671,533–1,153,002; RSE 4.18%–8.61%). After applying FHRF-R, the estimate decreased to IDR 1,055,405, and the RSE dropped to 15.50%, resulting in values more consistent with the regional pattern. This finding illustrates the model’s ability to correct extreme direct estimates arising from small sample sizes by borrowing strength across auxiliary variables and inter-area variability [1, 28].

tu_pian_32.png

The pattern illustrated by the case of Kota Manado is not isolated. Table 7 summarizes the ten districts/cities exhibiting the largest absolute differences between direct and SAE-FHRF-R point estimates. These areas are predominantly urban districts with very small agricultural household samples, ranging from 15 to 201 observations, which is far below the national average of 962 households per district/city. They are also characterized by high direct-estimate RSE values, indicating substantial sampling uncertainty.

Table 7. Districts with the largest changes in point estimates between direct estimation and SAE-FHRF-R

Note: SAE: Small area estimation; FHRF-C: Fay–Herriot with Random Forest correction; FHRF-R: Fay–Herriot with Random Forest replacement.

In these domains, the SAE-FHRF-R model produces systematic adjustments to point estimates, accompanied by marked reductions in RSE—for example, the 76% reduction observed in Kota Surakarta. By contrast, districts with adequate sample sizes and low direct-estimate RSEs exhibit negligible changes. This selective behavior explains the high overall correlation between direct and SAE-FHRF-R estimates (r = 0.87), which arises from the stability of well-sampled areas combined with targeted corrections in a limited number of small-sample urban domains.

Taken together, Figures 4 and 5 along with the summary in Table 7 indicate that substantial differences between direct and FHRF-R estimates arise almost exclusively in small-sample urban districts (area codes with third digit “7”), where agricultural households are scarce. Several cities with RSEs above the 25% threshold—such as Surakarta, Manado, Bandung, and Cirebon—exhibit substantial reductions in RSE under FHRF-R. In contrast, cities with moderate direct RSEs display smaller or mixed changes. In contrast, the majority of districts with adequate sample sizes and low direct-estimate RSEs remain largely unchanged, resulting in only localized differences in the spatial maps (Figure 6).

tu_pian_33.png

(a) Direct estimates

tu_pian_34.png

(b) SAE-FHRF-R estimates

tu_pian_35.png

Figure 7. Random Forest variable importance measured by mean decrease in accuracy (top ten)

In addition to improving precision, FHRF-R yields interpretable information from its variable-importance results. As shown in Figure 7, the share of villages with internet cafés was the strongest predictor of agricultural household expenditure, suggesting the role of digital infrastructure in local economic conditions [39]. This information complements the model’s predictive advantages by providing empirically grounded insights that can support targeted policy interventions [14, 21].

Overall, the empirical results showed that SAE-FHRF-R outperforms both direct estimation and FH–EM under nonlinear conditions, yielding higher precision and stable performance in small-sample urban areas. These findings are consistent with prior literature [4, 14, 21] and are reinforced by recent developments in integrating machine learning methods into the SAE framework [4, 10]. The variable-importance results further identify the most influential determinants of agricultural household expenditure, enhancing the model’s policy relevance.

Consistent with this pattern, the spatial maps in Figure 6 exhibit only localized differences between direct estimates and FHRF-R estimates. This pattern occurs because most adjustments are concentrated in small urban domains with high RSEs, whereas the majority of districts—with already low RSEs—receive minimal correction.