Sub-national heterogeneity in the time-varying reproduction number during the 2026 Bundibugyo virus disease outbreak in the Democratic Republic of the Congo: a hierarchical Bayesian analysis

National-level estimates of the time-varying reproduction number (Rt) for the 2026 Bundibugyo virus disease (BDBV) outbreak in the Democratic Republic of the Congo (DRC) have converged on a value close to the epidemic threshold since early August 2026, consistent with independent joint Bayesian renewal-model estimates. A single national Rt, however, can obscure divergent sub-national epidemic trajectories, particularly across a five-province outbreak in which provinces range from a declining original epicentre to recently-seeded fronts. We estimated Rt at national, provincial, and, where case volume allowed, health-zone level, using both a standard sliding-window (Cori) estimator and a hierarchical Bayesian renewal model with partial pooling across spatial units, fitted by Hamiltonian Monte Carlo (No-U-Turn Sampler). Provincial estimates diverged materially from the national trend: as of the week of 6 August 2026, Ituri, the outbreak's original epicentre, had a hierarchical median Rt of 0.91 (95% credible interval [CrI] 0.68 - 1.24), while Nord-Kivu (1.23, [0.87 - 1.72]) and Haut-Uele (1.79, [1.24 - 2.67]) remained above threshold. Health-zone disaggregation, feasible only in Ituri and Nord-Kivu given case volume, showed this provincial picture itself masked further heterogeneity: in Ituri, the zone where the outbreak began (Mongbwalu) had clearly declined (Rt 0.36, [0.17 - 0.74]) while the two largest zones by cumulative case count (Bunia, Rwampara) remained at or above threshold; in Nord-Kivu, elevated transmission was concentrated in a single zone (Katwa, Rt 1.39, [0.85 - 1.99]) while a comparably-sized zone (Butembo) had already declined (0.64, [0.23 - 1.38]). An initial disagreement between the sliding-window and hierarchical provincial estimates was traced to a data-reconstruction artefact (forward-filling, rather than interpolating, multi-day gaps in health-zone reporting) rather than a genuine methods disagreement, and resolved once corrected. The hierarchical model's dispersion structure, calibration, and sensitivity to the generation-interval assumption were each checked explicitly; a shared (non-province-specific) dispersion parameter was retained on the basis of negligible predictive difference (PSIS-LOO), the model achieved 95.0% pooled 95% posterior-predictive interval coverage, and the province ranking was unchanged across a generation-interval sensitivity grid (Spearman {rho} = 1.0). Aggregation masks meaningful heterogeneity in transmission intensity at every spatial resolution examined; response prioritisation based on a single national or even provincial Rt risks directing attention away from the specific zones where transmission remains supercritical.