  Estimates of national immunization coverage are crucial for guiding policy
and decision-making in national immunization programs and setting the global
immunization agenda. WHO and UNICEF estimates of national immunization coverage
(WUENIC) are produced annually for various vaccine-dose combinations and all
WHO Member States using information from multiple data sources and a
deterministic computational logic approach. This approach, however, is
incapable of characterizing the uncertainties inherent in coverage measurement
and estimation. It also provides no statistically principled way of exploiting
and accounting for the interdependence in immunization coverage data collected
for multiple vaccines, countries and time points. Here, we develop Bayesian
hierarchical modeling approaches for producing accurate estimates of national
immunization coverage and their associated uncertainties. We propose and
explore two candidate models: a balanced data single likelihood (BDSL) model
and an irregular data multiple likelihood (IDML) model, both of which differ in
their handling of missing data and characterization of the uncertainties
associated with the multiple input data sources. We provide a simulation study
that demonstrates a high degree of accuracy of the estimates produced by the
proposed models, and which also shows that the IDML model is the better model.
We apply the methodology to produce coverage estimates for select vaccine-dose
combinations for the period 2000-2019. A contributed R package {\tt imcover}
implementing the No-U-Turn Sampler (NUTS) in the Stan programming language
enhances the utility and reproducibility of the methodology.
