0
$\begingroup$

I'm working with a bimodal data set and need to split this bimodal data to 1) estimate the proportion of data points in each distribution and 2) the mean of each distribution.

I have approached this a couple ways, 1) by fitting mixed models, and 2) by estimating a cutoff value between the peaks. Fitting mixed models using a few packages in R (mixtools, mclust, and mixsmsn) worked well for me for a long time. However, I have new data where all of these methods have failed to fit my data. I've also tried various log transformations of my data, but I still can't get accurate fits.

I'm stumped. I've tried experimenting with providing estimated parameters for the model fits in mixtools, mclust, and mixsmsn but still no luck fitting. I think that every method I have tried thinks that the lower bimodal peak in the trouble data is just noise, but I know it is not!

Histograms of bimodal data: enter image description here

link to bimodal data (no issue fitting)

link to trouble data

$\endgroup$

1 Answer 1

1
$\begingroup$

The assumption of those mixed models (at least mclust and mixtools::normalmixEM2comp) is that the data is a mixture of Gaussians. I think in your "trouble" case, the first peak is just not Gaussian enough. So, if you want to stay with approach 1, it might be possible to model your data as a mixture of a Gaussian and a point-probability. But this is not easy, I'm not sure R has ready-made functions for that purpose.

Looking at the histograms, it seems you have two nice peaks, well-separated with a big valley. In that case, it might be easier to look for the valley than to look for a mixture (your approach number 2). For example, here is a very simple and naive approach of smoothing the histogram and looking for a negative peak:

trouble_data <- read.csv("https://raw.githubusercontent.com/jenniferbethann/bimodaldata/main/troubledata.csv",
                header=TRUE)$trouble_data

troub <- log10(trouble_data)

hist(troub, freq = FALSE, breaks = 100)
smtroub <- density(troub)
lines(smtroub, col="darkblue")
abline(v = smtroub$x[which(peakPick::peakpick(-smtroub$y,neighlim = 1))],col="red")

Created on 2020-12-09 by the reprex package (v0.3.0)

Better approaches should exist, in particular those developed for image analysis. For example you can take a look at the methods available in ImageJ. The package autothresholdr provides them in R but requires the data to be integers and has problems if it doesn't start at 0. We can try to run it by transforming the data a bit more:

trouble_data <- read.csv("https://raw.githubusercontent.com/jenniferbethann/bimodaldata/main/troubledata.csv",
                         header=TRUE)$trouble_data

troub <- as.integer(1000*log10(trouble_data))
troub <- troub - min(troub)

hist(troub, freq = FALSE, breaks = 100)

# my previous method
smtroub <- density(troub)
lines(smtroub, col="darkblue")
abline(v = smtroub$x[which(peakPick::peakpick(-smtroub$y,neighlim = 1))],col="red")


methods <- c("IJDefault", "Huang", "Huang2", "Intermodes", "IsoData", "Li", 
             "MaxEntropy", "Mean", "MinErrorI", "Minimum",
             "Moments", "Otsu", "Percentile", "RenyiEntropy", "Shanbhag",
             "Triangle", "Yen")
auto_thresholds <- sapply(methods,
                          function(x) tryCatch(autothresholdr::auto_thresh(troub, x), 
                                               error=function(e) NA))

abline(v=auto_thresholds,
       col="green")

Created on 2020-12-09 by the reprex package (v0.3.0)

So it looks like the InterModes method does perform well, the others not so much.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.