On cross-validation for small area estimators

Let $U_{i}\subset\{1,2,...,N\}$ denote the set of indices for the units that belong to the $i$-th subpopulation or area. We focus on binary outcomes in this paper, though the method can be generalized to other types of variables too. Let $y_{j}\in\{0,1\}$ denote the response value from the $j$-th unit and $s[j]\in\{1,...,M\}$ denote which subpopulation the $j$-th unit belongs to. From a probabilistic survey, we obtain a sample $S\subset U_{1}\cup U_{2}\cup\cdots\cup U_{M}$ of the units, with their design weights $w_{j}=1/\pi_{j}$, where $\pi_{j}$ is the probability that the $j$-th unit is sampled. To simplify the notation when comparing different sample splitting strategies, we denote a generic partition of a sample $S$ into a training set $A$ and a validation set $B=S\setminus A$, and use $\hat{\theta}_{A}$ and $\hat{\theta}_{B}$ to denote estimators based on the two subsamples. Specifically for $K$-fold cross-validation, we use $\hat{\theta}^{(k)}$ and $\hat{\theta}^{(-k)}$ to denote the estimators based on the $k$-th held-out dataset and the remaining samples respectively. As a working example, we consider the stratified multi-stage design adopted by most DHS surveys, where the primary sampling units (PSUs) are enumeration areas (EAs), or clusters, and within each PSU, the secondary sampling units (SSUs) are households.

Let $\theta_{i}$ denote the true prevalence in area $i$ and $\bm{\theta}=(\theta_{1},...,\theta_{M})$. When data are sufficient for direct estimation, the Hájek estimator (Hájek, 1971) for subpopulation prevalence is given by

Area-level models then proceed by smoothing direct estimates $\hat{\bm{\theta}}^{w}$ and the associated design-based variance, whereas unit-level approaches specify a parametric distribution for $y_{j}$ directly. Our objective is to evaluate the resulting estimates $\hat{\bm{\theta}}$ regardless of the underlying model.

For most SAE applications, a natural metric to score the estimator is the weighted squared error loss, $L(\hat{\bm{\theta}},\bm{\theta})=\sum_{i=1}^{M}q_{i}(\hat{\theta}_{i}-\theta_{i})^{2},$ for some weights $q_{i}$. The choice of the aggregation weights is contextual. They may represent population proportions by area, or simply be set to $1/M$ to weigh all areas equally. Evaluating $\hat{\bm{\theta}}$ based on the loss function instead of model-induced likelihood enables model comparison across different model classes. Similar principles have been adopted by Kuh et al. (2024) in the context of multilevel regression and poststratification.

This squared error loss is closely related to the mean squared error (MSE), a topic extensively studied in the SAE literature. Specifically, our ultimate target of inference for model $\mathcal{M}$ is the weighted average of area-specific conditional MSEs,

where $\bm{\theta}$ is considered fixed and unknown, and the expectation is taken over the distribution of sampling units. The term $\mathbb{E}[(\hat{\theta}_{i}^{\mathcal{M}}-\theta_{i})^{2}]$ is usually referred to as the conditional MSE in the literature (Rivest & Belmonte, 2000; Lohr & Rao, 2009; Rao & Molina, 2015). The majority of the theoretical results on estimating MSE, however, focus on the unconditional case where the expectation is further taken over the distribution of $\bm{\theta}$, assuming area-level Gaussian random effects (Prasad & Rao, 1990; Datta & Lahiri, 2000). To the best of our knowledge, the score proposed in this paper is the first cross-validation-based estimator of the conditional MSE without imposing distributional assumptions on the random effects.

Following the multi-stage survey design, we randomly assign SSUs, i.e., households, within each cluster to one of the $K$ folds. This creates $K$ folds of data each containing all clusters in the survey but only $1/K$ of the households within each cluster. The survey design is preserved within each fold, and the selection probability for each unit only needs to be adjusted by $1/K$. Given a partition, the naive estimator of the squared error loss is

where $\hat{\theta}_{i}^{(-k)}$ denotes the model-based prediction for area $i$ obtained using all data except the $k$-th fold, and $\hat{\theta}_{i}^{(k),w}$ denotes the direct estimate for area $i$ computed from the $k$-th held-out fold. The partitioning step can also be repeated and averaged to stabilize the estimated score.

Stratified cross-validation under complex survey designs has been investigated in the survey literature, usually for the purpose of variable selection in regression models (Iparragirre et al., 2023; Wieczorek et al., 2022). We note that our goal here is different. Under $K$-fold splitting, the estimators are compared with the direct estimates of the held-out samples, instead of predicting for individual survey responses. This allows the resulting score to be comparable for both area- and unit-level models.

Importantly, validation approaches based on sample splitting target the conditional MSE at a reduced sample size, i.e., $\text{MSE}_{A}^{\mathcal{M}}=\sum_{i=1}^{M}q_{i}\mathbb{E}[(\hat{\theta}_{A,i}^{\mathcal{M}}-\theta_{i})^{2}]$, where $A$ is a sample of size $(K-1)/K$ of $S$ approximately, rather than $\text{MSE}^{\mathcal{M}}$ in Equation 2. To differentiate the two targets, we refer to $\text{MSE}^{\mathcal{M}}$ as the oracle full-sample MSE and $\text{MSE}_{A}^{\mathcal{M}}$ the oracle training MSE. A similar gap in the target estimand was also observed in Kawano et al. (2026) and it generally cannot be estimated. For the purpose of model comparison, our goal is to preserve model ranking under the original design. It is therefore desirable to keep the training sets close to the full data distribution. For multi-stage designs, one may choose to stratify at different stages to form the partitions. In the Supplementary Materials, we include one alternative partitioning strategy where clusters are assigned to the $K$ folds instead of households, and show that such approach performs similarly to cross-validation of SSUs when sample size is large, but is more biased when sample size is small, due to a larger gap between the training oracle and full-sample oracle estimand.

Given a sample splitting strategy, the naive CV score is still biased for $\text{MSE}_{A}$, because the held-out direct estimate carries sampling variation. Theorem 3.1 identifies the bias structure for sample splitting estimators in general.

While $\hat{m}_{S}$ provides an unbiased estimator of the training MSE, the conditional biases, $b_{i}$ and $d_{i}^{\mathcal{M}}$, are not identifiable without knowing the true $\bm{\theta}$. It is important to note that while the direct estimator is design unbiased, splitting a survey into two sets results in the realized error being carried over to all resulting sub-samples. Thus $b_{i}$ is nonzero given a finite sample $S$. This leads to an unidentifiable bias when evaluating the difference of MSEs across models. Theorem 3.3 provides an adjusted score and characterizes the remainder bias term for cross-validation estimators.

Among the two correction terms, $v_{i}$ is a constant across models and $c_{i}^{\mathcal{M}}$ is analogous to the covariance penalty in Efron (2004), reflecting the shared survey structures between the training and validation samples. Theorem 3.3 also shows that the residual bias in the score difference is driven by the covariance between the conditional bias of the direct estimate in the validation set, $b_{i}$, and the model difference in the training set, $d_{i}^{\mathcal{M}_{1}}-d_{i}^{\mathcal{M}_{2}}$. The bias is large when, across hypothetical replicate surveys, the sampling error in the direct estimator for area $i$ is systematically associated with the difference in model biases. This occurs, for instance, when one model tracks the direct estimates more closely than the other. When there is only a single sample $S$, the bias $e_{i}$ cannot be identified. Therefore, when two models have similar conditional MSE, the oracle MSE difference in Equation (12) is small and the remainder $e_{i}$ may dominate the score difference, making the comparison uninformative. To account for the unidentifiable bias, we establish a finite sample bound for $|e_{i}$ in Theorem 3.4, and show that $e_{i}\to 0$ as data accumulates in Theorem 3.6.

Theorem 3.4 gives rise to a plug-in estimator of the upper bound of the absolute bias. Under the stratified cross-validation scheme,

The first term $\widehat{\mbox{var}}(\hat{\theta}_{S,i}^{w})$ is the standard design-based variance estimator for the direct estimate based on the full sample. When considering Bayesian models, the posterior variance, $\mbox{var}_{\text{post}}(\theta_{S,i}^{\mathcal{M}})$ provides a conservative plug-in estimator for $\widehat{\mbox{var}}(\hat{\theta}_{S,i}^{\mathcal{M}})$, since it exceeds the sampling variance of the posterior mean by accounting for residual model uncertainty in addition to sampling variability.

For the estimation of the aggregated squared error, the aggregated adjusted CV score is $\widehat{\mbox{score}}_{q}^{\mathcal{M}}=\sum_{i=1}^{M}q_{i}\widehat{\mbox{score}}_{i}^{\mathcal{M}}$, with the aggregated bias $e_{q}=\sum_{i=1}^{M}q_{i}e_{i}$, satisfying

Thus $t_{q}=\sum_{i=1}^{M}q_{i}t_{i}$ serves as a threshold for deciding whether an observed aggregated score difference is large enough to support a ranking between $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$. The full procedure for $K$-fold cross-validation of small area estimators is summarized in Algorithm 1.

Finally, the preceding results provide finite-sample tools for model comparison. A natural question is whether the remainder $e_{i}$ vanishes as data accumulate. Theorem 3.6 shows that the adjusted score is asymptotically consistent for the conditional MSE difference under mild conditions.

To illustrate the pitfalls of sample splitting strategies that fail to maintain the original survey design, we examine a popular validation scheme in the literature using leave-one-area-out (LOAO) cross-validation. Under LOAO, the model is fitted on data from $M-1$ areas and then used to predict the prevalence in the held-out area. Let $\hat{\theta}_{i}^{(-i)}$ denote the predicted prevalence for the $i$-th area using the rest of the data, an analogous score evaluating the squared error loss is

LOAO validation has been informally used in the SAE literature (Steorts et al., 2020; Wakefield et al., 2025). For unit-level models, where individual observations are nested within areas, LOAO is a form of block cross-validation. For area-level models, LOAO corresponds directly to the standard leave-one-out cross-validation (LOO-CV) and can be evaluated using Pareto smoothed importance sampling without refitting the model $M$ times (Vehtari et al., 2016). The computational efficiency makes it a widely used diagnostic tool to assess Bayesian area-level models .

Although conceptually simple and straightforward to implement, LOAO is an extrapolation estimator targeting a different estimand, $\mathbb{E}[(\hat{\theta}_{i}^{(-i)}-\theta_{i})^{2}]$. The shifted target can lead to reversed rankings when comparing smoothing models in the typical SAE setting.

The proof of Proposition 3.7 follows from standard MSE decomposition. The gap between the two MSE estimand manifests differently across models. It is easy to see that for a strongly over-smoothed model, the gap is close to $0$ as the model estimates are not sensitive to dropping data from a single area. However, for a well calibrated model, both the estimation variance and bias of the estimator for the $i$-th area are expected to increase after removing all data from the $i$-th area, thus making the extrapolation MSE larger than the smoothing MSE. Therefore, LOAO systematically penalizes flexible models more than over-smoothed models. We numerically illustrate this bias in Section 4.2.

The 2024 Zambia DHS adopted a two-stage stratified cluster design with a sampling frame based on the 2022 Census of Population and Housing of the Republic of Zambia. The population comprises $36,770$ enumeration areas (EAs). Zambia has $10$ provinces, which are further divided into $115$ districts. In the first stage, each of the province was stratified by urban and rural residence, yielding $20$ strata in total. Within each stratum, EAs were selected using probability proportional to size (PPS) sampling, resulting in $545$ clusters in total. In the second stage, within each selected cluster, households were selected through equal-probability systematic sampling, with a fixed take of $25$ households per cluster, yielding a target sample of $13,625$ households. As a working example, we consider estimating the prevalence of female literacy in Zambia. We focus our analysis on prevalence estimation at the province level. The Supplementary Materials contain additional numerical results comparing models at the district level.

To illustrate the proposed cross-validation framework with models commonly used in practice, we consider three candidate estimators. At the area level, we consider Fay-Herriot models with

where $\hat{\phi}_{i}^{w}=\mbox{logit}(\hat{\theta}_{i}^{w})$ is the transformed direct estimates, and $\hat{V}_{i}$ is the associated asymptotic variance estimate of $\hat{\theta}_{i}^{w}$, and $u_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma_{u}^{2})$ are random effects.

At the unit level, Bernoulli likelihood for the binary survey responses is equivalent to a cluster-level binomial model. For a survey consisting of $C$ clusters, we let $c[j]$ be the cluster index of the $j$-th individual and $i[c]$ be the area index of the $c$-th cluster. Let $Y_{c}$ and $n_{c}$ denote the number of positive responses and the number of sampled units from the $c$-th cluster, respectively, i.e., $Y_{c}=\sum_{j:c[j]=c}y_{k}$ and $n_{c}=\sum_{j:c[j]=c}1$. We consider a cluster-level model with a beta-binomial likelihood (Dong & Wakefield, 2021; Dong et al., 2026; Wakefield et al., 2025),

In this formulation, $p_{c}$ is the latent prevalence for the $c$-th cluster, and $d$ is the overdispersion parameter that captures the additional within-cluster variation. The cluster-level prevalence is then linked to area-level prevalence via $\mbox{logit}(p_{c})=\alpha+u_{i[c]},$ with random effects $u_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma_{u}^{2})$.

For the rest of the section, we consider three candidate models that differ in their model class and prior specification. For all three models, we adopt the penalised complexity (PC) prior (Simpson et al., 2017), where $PC(U,\alpha)$ prior for the standard deviation $\sigma$ specifies that $Pr(\sigma>U)=\alpha$.

• •

  $\mathcal{M}_{1}$: Unit-level Beta-Binomial model with a default $PC(1,0.01)$ prior on $\sigma_{u}$.

• •

  $\mathcal{M}_{2}$: Area-level Fay-Herriot model with the same priors for random effects as in $\mathcal{M}_{1}$.

• •

  $\mathcal{M}_{3}$: Same model as in $\mathcal{M}_{1}$, but with a highly informative $PC(0.01,0.01)$ prior on $\sigma_{u}$.

All models are fitted using the R package surveyPrev (Dong et al., 2024), employing INLA (Lindgren & Rue, 2015; Rue et al., 2009) in posterior computation. Figure 1 compares the direct estimates and estimates from the candidate models. These estimates will be formally evaluated using the proposed method in Section 4.3.

We first consider a design-based simulation study that mimics the 2024 Zambia DHS. Since the 2024 Zambia DHS was based on the sampling frame from the 2022 census, we constructed a 300m $\times$ 300m population grid using the 2022 population surface from WorldPop (WorldPop, 2018). To construct the sampling frame, we ranked pixels by population within each province and retained the most populated pixels to match the provincial EA counts from the 2024 DHS sampling frame, yielding $36{,}770$ gridded pixels. We treat these pixels as the master sampling frame of our study. For simplicity, we ignore urban/rural status of the pixels in the simulation and consider synthetic surveys stratified by province.

While the analysis in this paper focuses on prevalence estimation at the province level, we construct the synthetic population in this simulation study using smoothed district-level estimates of women’s literacy rate, obtained by fitting a Fay-Herriot model on the 2024 Zambia DHS. This introduces realistic within-province heterogeneity in our simulated data. For a given cluster $c$ with population $N_{c}$ in a district with assumed true prevalence $p$, we generate event count $Y_{c}$ from a Beta-Binomial distribution with mean $N_{c}p$ and variance $N_{c}p(1-p)(1+(N_{c}-1)d)$, where $d=0.00001$ introduces additional overdispersion within the district. The simulated finite population prevalence in the $i$-th province is then $\theta_{i}=\sum_{c\in i}Y_{c}/\sum_{c\in i}N_{c}$.

Given the synthetic population, we draw 50 replication surveys using a two-stage stratified design, where $50$ clusters are sampled in each stratum and $30$ households are sampled in each cluster. Figures 2 displays, for a single replicate, point estimates and 90% uncertainty intervals produced by the candidate models. Direct estimation, $\mathcal{M}_{1}$, and $\mathcal{M}_{2}$ all produce estimates close to the population truth across all Admin-1 areas. By contrast, $\mathcal{M}_{3}$ shows excessive shrinkage, though with overlapping interval estimates compared to others.

We computed the adjusted cross-validation scores of all models using $5$-fold stratified cross-validation. Figure 3 summarizes the adjusted CV score comparisons of the three models across the $50$ synthetic surveys. When comparing $\mathcal{M}_{3}$ with $\mathcal{M}_{1}$, the adjusted CV score was able to consistently identify that $\mathcal{M}_{3}$ has larger MSE and the score differences all exceed the error bound, $t_{q}$, providing decisive evidence in favor of $\mathcal{M}_{1}$. The comparison between $\mathcal{M}_{2}$ and $\mathcal{M}_{1}$ is more difficult as the oracle differences are much smaller in scale. The adjusted CV score differences mostly have the same sign as the oracle difference, but in all replicates, the estimated score differences were between $\pm t_{q}$, thus suggesting the comparison between $\mathcal{M}_{2}$ and $\mathcal{M}_{1}$ is inconclusive, as the available information in the data does not overcome the unestimable bias components of the conditional MSE.

In this simulation setting, the relative performance of the candidate models remains unchanged between the full-sample and oracle training MSE. When sample sizes are smaller, the two estimands may lead to different model rankings, which we investigate in the Supplementary Materials. We also compare different sample-splitting strategies and the choice of $K$ in the Supplementary Materials.

To gain additional insight on the model performance, we also examine the area-specific scores in Figure 4. For the comparison between $\mathcal{M}_{3}$ and $\mathcal{M}_{1}$, the over-smoothed $\mathcal{M}_{3}$ departs noticeably from the other models in Copperbelt, Luapula, Lusaka, and Northern provinces, producing absolute score differences that are well above the error bound. For comparing $\mathcal{M}_{2}$ with $\mathcal{M}_{1}$, area-specific score differences still remain tightly clustered near zero, with all points falling within the inconclusive region, except for Lusaka in a few occasions.

In summary, these results demonstrate that the CV score reliably ranks the models for the performance of subnational level prevalence estimation, while the error bound $t_{q}$ provides the principled threshold for separating genuine differences from noise-dominated ones. Additional diagnostics and visualizations are provided in Supplementary Materials.

We also compare the LOAO strategy in Figure 5 for $\mathcal{M}1$ and $\mathcal{M}3$ under leave-one-area-out (LOAO) validation. As discussed in Section 3.4, LOAO penalizes the over-smoothed model less. Indeed in this comparison, the procedure reached the opposite conclusion and ranked $\mathcal{M}3$ as the preferred model in every replicate, illustrating the bias described in Section 3.4.

We now return to the 2024 Zambia DHS survey with the indicator given by the percentage of women who are literate. The stratified cross-validation was conducted with the $20$ strata defined by province crossed with urban/rural residence. Results from 5-fold cross-validation is reported in Figure 6. Estimates from $\mathcal{M}_{3}$ is heavily shrunk to the overall mean. Between $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, we observe slightly more shrinkage from $\mathcal{M}_{1}$, and the shrinkage is more obvious in Lusaka and Copperbelt. The cross-validation procedure again identifies that $\mathcal{M}_{3}$ has larger aggregated score that exceeds the estimated threshold, compared with $\mathcal{M}_{2}$. The difference between $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, on the other hand, is a lot smaller and within the inclusive region.

A close examination of the area-specific scores in Figure 6 provides additional insight into the behavior of the models. When comparing $\mathcal{M}_{2}$ and $\mathcal{M}_{3}$, the score differences are negative and exceed the area-specific threshold in magnitude across all ten provinces, yielding a uniform preference for $\mathcal{M}_{2}$ over $\mathcal{M}_{3}$. When comparing $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, while the overall score differences are smaller than the potential bias in magnitude, there are two provinces, Lusaka and Copperbelt, with score differences exceeding their area-specific thresholds. This is consistent with what we observe with the estimated prevalence difference between the two models, and can lead to $M_{2}$ being preferred under another set of aggregation weights different from the population proportions in this analysis.