Bayesian multinomial small area estimation with post-stratification and stochastic benchmarking

Reliable estimation of categorical outcomes at detailed geographical levels is challenging when local survey samples are small or outcome categories are sparse. This study develops a Bayesian unit-level small area estimation framework for multinomial outcomes, integrating hierarchical multinomial logistic regression, multilevel regression and post-stratification (MRP), Iterative Proportional Fitting (IPF), survey-design-based bootstrap uncertainty propagation, and stochastic benchmarking to uncertain aggregate estimates. The framework accounts for uncertainty arising from model estimation, reconstruction of joint auxiliary-variable distributions, and benchmarking. The methodology is applied to depression severity across 54 Greek prefectures using the 2019 Greek Health Survey and auxiliary population information. The resulting area-level estimates were stable and internally coherent and generally more precise than direct survey estimates, particularly for less prevalent categories and areas with limited sample information. A design-aware Monte Carlo simulation with 200 replications showed generally small parameter bias, satisfactory recovery of area-level heterogeneity, and empirical coverage reasonably close to nominal levels. Benchmarking reduced mean squared error relative to unbenchmarked MRP, while deterministic and stochastic benchmarking yielded virtually identical point-estimation performance. Stochastic benchmarking additionally propagated uncertainty in the estimated aggregate benchmark, with slightly greater posterior uncertainty and modestly improved coverage. Overall, the empirical and simulation results support the finite-sample performance and practical applicability of the proposed multinomial small area estimation (SAE) framework under the conditions examined.