Molecular aggregation equilibria are described using finite lattice and mean field theoretical modeling strategies, both built upon a random mixture reference system. The resulting predictions are compared with each other for systems in which each aggregate consists of a central solute molecule whose first coordination shell can accommodate multiple bound ligands. Solute ligand (direct) and ligand ligand (cooperative) interactions are found to influence aggregate size distributions in qualitatively different ways, as direct interactions produce a shape-invariant transformation of the aggregate size distribution, whereas cooperative interactions can lead to a vapor liquid-like transformation. When half the ligand binding sites are filled, the corresponding aggregate size distributions are invariably unimodal in the absence of cooperative interactions, but when the latter interactions are attractive, the distributions are predicted to be bimodal below and unimodal above a critical temperature. Mean field and finite lattice predictions are found to be in globally good agreement with each other, except under near-critical conditions, and even there, the predicted average aggregate sizes and equilibrium constants are remarkably similar. Potential applications of these theoretical predictions to the analysis of experimental and molecular dynamics aggregation results are discussed.