# How can I estimate the probability that my observed data come from a bimodal population?

I have a data set representing the abundance of a protein in a population of cells. Based on our understanding of the biology behind this, I expect there to be two subpopulations - one in which this protein is abundant and one in which it is not. The observed data seem to fit this hypothesis well, and certainly look bimodal (I'll try to attach an example below). If we assume that any subpopulations are themselves normally-distributed, is there a way to quantify the probability that this data is multimodal as opposed to being unimodal?

Looking through similar questions elsewhere on this exchange, a mixture model seems to be designed to do this, but I have no idea how to implement that. Is this something that could realistically be accomplished by a statistics novice?

Many thanks, O

If you're assuming each sub-population is normally distributed (even in the case where there's only one), then you could use a likelihood-ratio test. Use a single Gaussian for the null hypothesis and the sum of two Gaussian random variables for the alternative hypothesis. Since these are nested models, Wilks' theorem says that twice the log likelihood ratio will be distributed approximately as a $\chi^2$ distribution in the limit, but if you don't have enough data, you can approximate this by Monte Carlo.
• In terms of the OP, Gaussian mixture modeling is implemented in the R package mclust, which in this case could be done on the log(Mean Intensity) values. That's fairly straightforward, but is based on the Gaussian assumption. The flexmix package provides for more general mixture modeling; in this case, one could model log(Total Intensity) as a function of Cell Area and look for 2 clusters with different slopes of that relation. – EdM Feb 4 '17 at 17:55