The data obtained by the recent modern sky surveys enable detailed studies of the stellar distribution in the multi-dimensional space spanned by spatial coordinates, velocity and metallicity, from the solar neighborhood all the way out to the outer Milky Way halo. While these results represent exciting observational breakthroughs, their interpretation is not simple. For example, traditional decomposition of the thin and thick disks predicts a strong correlation in metallicity and kinematics at $sim$1 kpc from the Galactic plane; however, recent SDSS--based work has demonstrated an absence of this correlation for disk stars. Instead, the variation of the metallicity and rotational velocity distributions can be modeled using non--Gaussian functions that retain their shapes and only shift as the distance from the mid--plane increases. To fully contextualize these recent observational results, a detailed comparison with sophisticated numerical models is necessary. Modern simulations have sufficient resolution and physical detail to study the formation of stellar disks and spheroids over a large baseline of masses and cosmic ages. We discuss preliminary comparisons of various observed maps and N--body model predictions and find them encouraging. In particular, the N--body disk models of Rov{s}kar et al. cite{Roskar 2008} reproduce a change of disk scale height reminiscent of thin/thick disk decomposition, as well as metallicity and rotational velocity gradients, while not inducing a correlation of the latter two quantities, in qualitative agreement with SDSS observations.