We present a multilevel Monte Carlo simulation method for analyzing multi-scale physical systems via a hierarchy of coarse-grained representations, to obtain numerically exact results, at the most detailed level. We apply the method to a mixture of size-asymmetric hard spheres, in the grand canonical ensemble. A three-level version of the method is compared with a previously studied two-level version. The extra level interpolates between the full mixture and a coarse-grained description where only the large particles are present-this is achieved by restricting the small particles to regions close to the large ones. The three-level method improves the performance of the estimator, at fixed computational cost. We analyze the asymptotic variance of the estimator and discuss the mechanisms for the improved performance.