We model the evolution of active galactic nuclei (AGN) by constructing their radio LFs. We used a set of surveys of varying area and depth, namely, the deep COSMOS survey of 1916 AGN sources; the wide, shallow 3CRR, 7C, and 6CE surveys, together containing 356 AGN; and the intermediate XXL-North and South fields consisting of 899 and 1484 sources, respectively. We also used the CENSORS, BRL, Wall & Peacock, and Config surveys, respectively consisting of 150, 178, 233, and 230 sources. Together, these surveys account for 5446 AGN sources and constrained the LFs at high redshift and over a wide range of luminosities (up to z~=3 and log(L/W/Hz){in}[22,29]). We concentrated on parametric methods within the Bayesian framework, which allowed us to perform model selection between a set of different models. By comparing the marginalised likelihoods and both the Akaike information criterion and the Bayesian information criterion, we show that the luminosity-dependent density evolution (LDDE) model fits the data best, with evidence ratios varying from 'strong' (>10) to 'decisive' (>100), according to the Jeffreys' interpretation. The best-fitting model gives insight into the physical picture of AGN evolution, where AGN evolve differently as a function of their radio luminosity. We determined the number density, luminosity density, and kinetic luminosity density as a function of redshift, and we observed a flattening of these functions at higher redshifts, which is not present in simpler models. We explain these trends by our use of the LDDE model. Finally, we divided our sample into subsets according to the stellar mass of the host galaxies in order to investigate a possible bimodality in evolution. We found a difference in LF shape and evolution between these subsets. All together, these findings point to a physical picture where the evolution and density of AGN cannot be explained well by simple models but require more complex models either via AGN sub-populations, where the total AGN sample is divided into sub-samples according to various properties, such as optical properties and stellar mass, or via luminosity-dependent functions.