A model for the simulation of ensembles of laser-driven Rydberg-Rydberg interacting multilevel atoms is discussed. Our hybrid approach combines an exact two-body treatment of nearby atom pairs with an effective approximate treatment for spatially separated pairs. We propose an optimized evolution equation based only on the system steady state, and a time-independent Monte Carlo technique is used to efficiently determine this steady state. The hybrid model predicts features in the pair-correlation function arising from multiatom processes which existing models can only partially reproduce. Our interpretation of these features shows that higher-order correlations are relevant already at low densities. Finally, we analyze the performance of our model in the high-density case.