Global species richness is a key biodiversity metric. Concerns continue to grow over its decline due to overexploitation and habitat destruction by humans. Despite recent efforts to estimate global species richness, the resulting estimates have been highly uncertain and often logically inconsistent. Estimates lower down either the taxonomic or geographic hierarchies are often larger than those above. Further, these estimates have been typically represented in a wide variety of forms, including intervals (a, b), point estimates with no uncertainty, and point estimates with either symmetrical or asymmetrical bounds, making it difficult to combine information across different studies. Here, we develop a Bayesian hierarchical approach to estimate global species richness (we estimate 22.02 m species; 95% highest posterior density (HPD) interval (10.43 million, 35.28 million)) that combines 50 estimates from published studies. The data mix of intervals and point estimates are reconciled using techniques from symbolic data analysis. This approach allows us to recover interval estimates at each species level, even when data are partially or wholly unobserved, while respecting logical constraints, and to determine the effects of estimation on the whole hierarchy of obtaining future estimates for particular taxa at various levels in the hierarchy.