All biological systems exhibit cell-to-cell variability. Frameworks exist for understanding how stochastic fluctuations and transient differences in cell state contribute to experimentally observable variations in cellular responses. However, current methods do not allow identification of the sources of variability between and within stable subpopulations of cells. We present a data-driven modeling framework for the analysis of populations comprising heterogeneous subpopulations. Our approach combines mixture modeling with frameworks for distribution approximation, facilitating the integration of multiple single-cell datasets and the detection of causal differences between and within subpopulations. The computational efficiency of our framework allows hundreds of competing hypotheses to be compared. We initially validate our method using simulated data with an understood ground truth, then we analyze data collected using quantitative single-cell microscopy of cultured sensory neurons involved in pain initiation. This approach allows us to quantify the relative contribution of neuronal subpopulations, culture conditions, and expression levels of signaling proteins to the observed cell-to-cell variability in NGF/TrkA-initiated Erk1/2 signaling.