Spatial aggregation with respect to a population distribution involves estimating aggregate population quantities based on observations from individuals. In this context, a geostatistical workflow must account for three major sources of : aggregation weights, fine scale variation, and finite population variation. However, these sources of aggregation error are commonly ignored, and the population instead treated as a fixed population density surface. We improve common practice by introducing a allowing aggregation models to account for aggregation error simply and transparently. This preserves aggregate point estimates while increasing their uncertainties. We compare the proposed and the traditional approach using two simulation studies mimicking neonatal mortality rate (NMR) data from the 2014 Kenya Demographic and Health Survey. In the traditional approach, undercoverage/overcoverage of interval estimates depends arbitrarily on the aggregation grid resolution, while the new approach is resolution robust. Differences between the aggregation approaches increase as an area's population decreases, and are particularly large at the second administrative level and finer, but also at the first administrative level for some population quantities. These findings are consistent with those of an application to the true NMR data. We demonstrate in a sensitivity analysis that burden estimates and their uncertainties are not robust to changes in population density and census information, while prevalence estimates and uncertainties seem stable.
