Authors: Katarzyna Reluga, Nicola Salvati, Setareh Ranjbar, Nikos Tzavidis, Mark van der Laan
Organizations: Chair of AI in Business and Economics, Humboldt University of Berlin, Germany · Department of Economics and Management, University of Pisa, Italy · Department of Psychiatry & Lausanne University Hospital, University of Lausanne, Switzerland · Department of Social Statistics and Demography & Southampton Statistical Sciences Research Institute, University of Southampton, UK · School of Public Health, University of California, Berkeley, USA
Small area estimation (SAE) traditionally borrows strength across domains through a prespecified working model. We introduce small area learning (SAL), a learner-agnostic framework that formulates SAE as a multi-target semi-supervised learning problem and treats outcome prediction as a nuisance task. Areas may differ through unrestricted covariate shift and structured conditional shift, with information shared through common and area-specific components. We develop standard and cross-fitted estimators for low-complexity and flexible learners, unifying synthetic, model-based, and fully residual-corrected estimation. For area means, we establish consistency and asymptotic normality and provide variance estimators and confidence intervals. We extend the theory to settings with only aggregated population covariates and characterize the resulting aggregation bias. SAL reveals a bias--variance trade-off: partial residual correction can reduce mean squared error, whereas full correction yields estimators robust to working-model misspecification, potentially at higher variance. Simulations yield RMSE reductions of 34%--49% relative to direct estimation, with the best configuration varying across settings. An empirical application favours a cross-fitted estimator using a flexible learner under the chosen validation criterion. The numerical studies also illustrate how shrinkage can produce overly optimistic confidence intervals at the expense of strong modeling assumptions.
Figures & tables
cell
scheme
learner
estimator
abs. bias
RMSE
gap (%)
coverage
S2_D25_n10
cf_p
bart_flexible
sae_naive
2.24
35.34
0.0
0.020
S2_D25_n10
cf_s
bart_flexible
sae_naive
2.23
35.39
0.2
0.022
S2_D25_n10
cf_p
bart
sae_naive
2.69
37.90
7.2
0.023
S2_D25_n10
sal
spline
sal
0.99
37.97
7.5
0.920
S2_D25_n10
—
—
direct
2.32
62.22
76.1
0.941
S2_D25_n20
cf_p
bart_flexible
sae_naive
0.83
24.92
0.0
0.038
Table 1: Scenario S2: the 3 combinations of learner, fitting scheme and estimator with the smallest RMSE for yˉd in each design cell, followed by the best fully corrected member where none ranks among them, and by the direct estimator. Absolute bias and RMSE are multiplied by 100; the gap is the percentage excess in RMSE over the best row of the cell. Coverage is of nominal 95% intervals with the tnd−1 quantile; for sae_naive it refers to the interval centred at the uncorrected predictor μ^d(fd) with half-width tnd−1,0.975fdσ^dnd−1/2 , which is not covered by Corollary 4.3 .
cell
scheme
learner
estimator
abs. bias
RMSE
gap (%)
coverage
S5_D25_n10
cf_p
bart_flexible
sae_naive
2.70
62.16
0.0
0.020
S5_D25_n10
cf_s
bart_flexible
sae_naive
2.69
62.23
0.1
0.023
S5_D25_n10
cf_s
bart
sae_naive
3.37
63.86
2.7
0.022
S5_D25_n10
sal
spline
sal
2.75
66.75
7.4
0.915
S5_D25_n10
—
—
direct
2.97
93.71
50.8
0.949
S5_D25_n20
cf_p
bart_flexible
sae_naive
2.56
43.08
0.0
0.039
Table 2: Scenario S5: the 3 combinations of learner, fitting scheme and estimator with the smallest RMSE for yˉd in each design cell, followed by the best fully corrected member where none ranks among them, and by the direct estimator. Absolute bias and RMSE are multiplied by 100; the gap is the percentage excess in RMSE over the best row of the cell. Coverage is of nominal 95% intervals with the tnd−1 quantile; for sae_naive it refers to the interval centred at the uncorrected predictor μ^d(fd) with half-width tnd−1,0.975fdσ^dnd−1/2 , which is not covered by Corollary 4.3 .
learner
estimator
scheme
mean
max
interval
rel. (%)
vs d. (%)
random_forest
sae
cf_s
363.805
578.517
no
1.95
52.9
random_forest
sae_cor
cf_s
422.121
694.576
yes
2.26
45.3
random_forest
sae
cf_p
375.677
599.741
no
2.01
51.3
random_forest
sae_cor
cf_p
400.450
694.210
yes
2.15
48.1
gradient_boosting
cf_a_bc
cf_p
377.511
921.081
yes
2.02
51.1
gradient_boosting
cf
cf_p
377.511
921.081
yes
2.02
51.1
Table 3: Best learner within each combination of estimator and fitting scheme, by the mean absolute distance ( 38 ) to the regional Horvitz–Thompson estimator over the 6 regions, with the maximum distance. Interval: whether the procedure supplies the interval of Corollary 4.3 . Rel.: mean distance as a percentage of the mean regional estimate. vs direct: percentage reduction in mean distance relative to the direct estimator. sae_cor is the bias-corrected version of the model-based predictor sae , which coincides with the fully corrected member of the same learner and fitting scheme. Below the rule: lmm , sae is the EBLUP under the nested-error model; ∗ its interval uses the parametric bootstrap MSE under that model.
learner
estimator
scheme
mean
max
interval
rel. (%)
vs d. (%)
random_forest
sae
cf_s
363.805
578.516
no
1.95
52.9
random_forest
sae_cor
cf_s
422.121
694.576
yes
2.26
45.3
random_forest
sae
cf_p
375.677
599.741
no
2.01
51.3
gradient_boosting
cf_a_bc
cf_p
377.511
921.081
yes
2.02
51.1
lmm
sae
sal
562.288
1031.664
yes ∗
3.01
27.1
—
direct
—
771.780
1267.491
yes
4.14
0.0
Table 4: The 3 best combinations of estimator and fitting scheme, each with its best learner, by the mean absolute distance ( 38 ) to the regional Horvitz–Thompson estimator over the 6 regions, with the maximum distance. Interval: whether the procedure supplies the interval of Corollary 4.3 . Rel.: mean distance as a percentage of the mean regional estimate. vs direct: percentage reduction in mean distance relative to the direct estimator. sae_cor is the bias-corrected version of the model-based predictor sae , which coincides with the fully corrected member of the same learner and fitting scheme. Below the rule: lmm , sae is the EBLUP under the nested-error model; ∗ its interval uses the parametric bootstrap MSE under that model.
rank
learner
estimator
scheme
mean estimate
mean CI length
1
random_forest
sae
cf_s
18540.459
—
—
random_forest
sae_cor
cf_s
18265.395
4895.149
2
gradient_boosting
sae
cf_s
18247.116
—
3
random_forest
sae
cf_p
18538.196
—
4
gradient_boosting
cf
cf_p
18060.852
5158.087
5
gradient_boosting
cf_a_bc
cf_p
18060.852
5158.087
Table 5: The 5 combinations of learner, estimator and fitting scheme closest to the regional Horvitz–Thompson estimator, ranked by the mean absolute distance ( 38 ) over the 6 regions, with the mean of the 41 provincial estimates and the mean length of the 95% intervals. The uncorrected model-based predictor sae has no interval of Corollary 4.3 ; sae_cor is its bias-corrected version, which coincides with the fully corrected member of the same learner and fitting scheme. Below the rule, for reference: lmm , sae , the EBLUP under the nested-error model ( ∗ interval from the parametric bootstrap MSE under that model), and the direct estimator.
Fig 1: Estimated mean equivalised household disposable income by province, ordered by sample size. Triangles: direct estimates ( 7 ) with design-based 95% intervals. Circles: the procedure closest to the regional target, the uncorrected predictor μ^d(fd) with a cross-fitted random forest (no interval). Squares: its bias-corrected version, the fully corrected cross-fitted random forest, with the intervals of Corollary 4.3 . Diamonds: the fully corrected cross-fitted gradient boosting learner with pooled folds, with the same intervals. Crosses: the EBLUP under the nested-error model, with intervals from its parametric bootstrap MSE.
S1
S2
S3
S4
S5
shared function m~0
( 42 )
( 43 )
( 44 )
( 44 )
( 44 )
deviation type
location
functional
functional
functional
functional
coefficient ratios s
—
(0.4,0.3,0.3)
(0,1,1)
(0,1,1)
(0,1,1)
covariate shift ψ
1
1
1.5
1.5
1.5
dependence τ
0
1
1.5
1.5
1.5
heterogeneity r
0.20
0.30
0.50
0.30
0.70
Table 6: Parameter values defining the five scenarios. Covariate dimensions are p=10 and q=2 throughout.
cell
D
Nd
nd
N
n
D25_n10
25
1000
10
25 000
250
D25_n20
25
1000
20
25 000
500
D50_n10
50
1000
10
50 000
500
D50_n20
50
1000
20
50 000
1000
Table 7: The four design cells, each run under all five scenarios. Sampling fractions nd/Nd are one or two per cent.
cell
scheme
learner
estimator
abs. bias
RMSE
gap (%)
coverage
S1_D25_n10
sal
lmm
sae_naive
1.05
28.33
0.0
0.020
S1_D25_n10
cf_strat
lmm
sae_naive
1.08
28.55
0.8
0.021
S1_D25_n10
cf_pooled
lmm
sae_naive
1.08
28.60
0.9
0.022
S1_D25_n10
sal
boosted_trees
sal
1.12
31.07
9.7
0.943
S1_D25_n10
—
—
direct
1.48
42.88
51.4
0.952
S1_D25_n20
sal
lmm
sae_naive
0.77
21.01
0.0
0.034
Table 8: Scenario S1: the 3 combinations of learner, fitting scheme and estimator with the smallest RMSE for yˉd in each design cell, followed by the best fully corrected member where none ranks among them, and by the direct estimator. Absolute bias and RMSE are multiplied by 100; the gap is the percentage excess in RMSE over the best row of the cell. Coverage is of nominal 95% intervals with the tnd−1 quantile; for sae_naive it refers to the interval centred at the uncorrected predictor μ^d(fd) with half-width tnd−1,0.975fdσ^dnd−1/2 , which is not covered by Corollary 4.3 .
cell
scheme
learner
estimator
abs. bias
RMSE
gap (%)
coverage
S3_D25_n10
cf_strat
bart
sae_naive
0.99
16.30
0.0
0.023
S3_D25_n10
cf_pooled
bart
sae_naive
1.00
16.34
0.3
0.023
S3_D25_n10
cf_pooled
bart_flexible
sae_naive
0.83
16.40
0.6
0.023
S3_D25_n10
sal
spline
sal
0.80
17.21
5.6
0.917
S3_D25_n10
—
—
direct
0.90
26.59
63.1
0.948
S3_D25_n20
cf_strat
bart_flexible
sae_naive
0.70
11.25
0.0
0.040
Table 9: Scenario S3: the 3 combinations of learner, fitting scheme and estimator with the smallest RMSE for yˉd in each design cell, followed by the best fully corrected member where none ranks among them, and by the direct estimator. Absolute bias and RMSE are multiplied by 100; the gap is the percentage excess in RMSE over the best row of the cell. Coverage is of nominal 95% intervals with the tnd−1 quantile; for sae_naive it refers to the interval centred at the uncorrected predictor μ^d(fd) with half-width tnd−1,0.975fdσ^dnd−1/2 , which is not covered by Corollary 4.3 .
cell
scheme
learner
estimator
abs. bias
RMSE
gap (%)
coverage
S4_D25_n10
cf_strat
spline
sae_naive
2.98
36.43
0.0
0.021
S4_D25_n10
sal
spline
sal
2.68
36.45
0.1
0.912
S4_D25_n10
sal
spline
sal_a_bc
2.68
36.45
0.1
0.912
S4_D25_n10
—
—
direct
2.67
77.62
113.1
0.943
S4_D25_n20
cf_strat
bart_flexible
sae_naive
1.96
22.14
0.0
0.042
S4_D25_n20
cf_pooled
bart_flexible
sae_naive
2.00
22.24
0.4
0.043
Table 10: Scenario S4: the 3 combinations of learner, fitting scheme and estimator with the smallest RMSE for yˉd in each design cell, followed by the best fully corrected member where none ranks among them, and by the direct estimator. Absolute bias and RMSE are multiplied by 100; the gap is the percentage excess in RMSE over the best row of the cell. Coverage is of nominal 95% intervals with the tnd−1 quantile; for sae_naive it refers to the interval centred at the uncorrected predictor μ^d(fd) with half-width tnd−1,0.975fdσ^dnd−1/2 , which is not covered by Corollary 4.3 .
Evaluating an AI system requires disaggregated assessment, as performance varies across domains such as benchmark task types or conversation types in deployed agents. Exhaustive testing is expensive, so evaluation rests on a sample of labeled units. We treat the evaluation set as a finite population and seek accurate point and interval estimates of each domain mean. Direct estimators, including prediction-powered inference (PPI), use only a domain's own labels and are imprecise where labels are few. Small area estimation addresses this problem, and we build on it to develop an integrated workflow for estimation and validation. For estimation, we propose prediction-powered smoothing (PP-S), a Bayesian model fit to each domain's prediction-powered estimate, with an extension that borrows strength across a reporting taxonomy (PP-TS). For validation, we derive a new, approximately unbiased design-based cross-validation score for choosing among direct and smoothed estimators. We study a curated benchmark with verifiable grading and deployed agent traffic graded by humans, each with every outcome observed. In both, the proposed estimators improve on the direct estimators in point and interval estimation, with near-nominal coverage. At the same sampling budget, our score selects as well as an independent validation sample does and estimates the selected estimator's error far more accurately.
Sho Kawano, Zehang Richard Li, Paul A. Parker
Department of Statistics, University of California, Santa Cruz, CA, USA
Estimating how much an intervention helps a given individual the conditional average treatment effect (CATE) is increasingly central to decision-making in medicine, economics, and policy, where an estimate is most useful when accompanied by a calibrated uncertainty interval. We study the few-placebo regime, in which one treatment arm is much smaller than the other, as arises in unequal-allocation trials and small-holdout A/B tests. The standard estimator in this setting is the X-Learner, and a natural way to obtain credible intervals is to make its second stage Bayesian. We show that these intervals under-cover: they contain the true effect less often than their nominal level. We trace this to a structural cause the X-Learner's regression target inherits the bias of a nuisance model fitted to the small arm, so the posterior is centered away from the true effect and we find that the standard remedy, regressing an orthogonal doubly-robust score, is also unreliable here, since the regime's limited overlap leaves the estimator either highly variable or, once stabilized, biased once more. Both consequences reflect a pattern that extends beyond causal inference: a separately estimated variance is attached to a point estimate of a hard-to-learn quantity, and the point estimate's bias is not captured by that variance. We propose GP-CATE, which models each arm's outcome surface with a Gaussian process, so the scarce arm's uncertainty enters the posterior directly rather than as an unmodelled bias. Across synthetic and semi-synthetic benchmarks, GP-CATE attains calibrated coverage where the estimators we compare against including Causal Forest and BART do not, at the cost of intervals that are appropriately wide when the data are uninformative.
When outcome data are expensive or onerous to collect, scientists increasingly substitute predictions from machine learning and AI models for unlabeled cases, a process which has consequences for downstream statistical inference. While recent methods provide valid uncertainty quantification under independent sampling, real-world applications involve missing at random (MAR) labeling and spatial dependence. For inference in this setting, we propose a doubly robust estimator with cross-fit nuisances. We show that cross-fitting induces fold-level correlation that distorts spatial variance estimators, producing unstable or overly conservative confidence intervals. To address this, we propose a jackknife spatial heteroscedasticity and autocorrelation consistent (HAC) variance correction that separates spatial dependence from fold-induced noise. Under standard identification and dependence conditions, the resulting intervals are asymptotically valid. Simulations and benchmark datasets show substantial improvement in finite-sample calibration, particularly under MAR labeling and clustered sampling.
Stephen Salerno, Zhenke Wu, Tyler McCormick
Fred Hutch Cancer Center · University of Michigan · University of Washington