Many statistical models have been proposed to analyse small area disease data with the aim of describing spatial variation in disease risk. In this paper, we propose a Bayesian hierarchical model that simultaneously allows for risk estimation and cluster identification. Our model formulation assumes that there is an unknown number of risk classes and small areas are assigned to a risk class by means of independent allocation variables. Therefore, areas within each cluster are assumed to share a common risk but they may be geographically separated. The posterior distribution of the parameter representing the number of risk classes is estimated using a novel procedure that combines its prior distribution with an efficient estimate of the marginal likelihood of the data given this parameter. An extension of the model incorporating covariates is also shown. These covariates may incorporate additional information on the problem or they may account for spatial correlation in the data. We illustrate the performance of the proposed model through both a simulation study and a case study of reported cases of varicella in the city of Valencia, Spain.