How to test if my distribution is multimodal?

Question

When I plot a histogram of my data, it has two peaks:

Does that mean a potential multi-modal distribution? I ran the dip.test in R (library(diptest)), and the output is:

D = 0.0275, p-value = 0.7913

I can conclude that my data have a multi-modal distribution?

DATA

10346 13698 13894 19854 28066 26620 27066 16658  9221 13578 11483 10390 11126 13487 
15851 16116 24102 30892 25081 14067 10433 15591  8639 10345 10639 15796 14507 21289 
25444 26149 23612 19671 12447 13535 10667 11255  8442 11546 15958 21058 28088 23827 
30707 19653 12791 13463 11465 12326 12277 12769 18341 19140 24590 28277 22694 15489 
11070 11002 11579  9834  9364 15128 15147 18499 25134 32116 24475 21952 10272 15404 
13079 10633 10761 13714 16073 23335 29822 26800 31489 19780 12238 15318  9646 11786 
10906 13056 17599 22524 25057 28809 27880 19912 12319 18240 11934 10290 11304 16092 
15911 24671 31081 27716 25388 22665 10603 14409 10736  9651 12533 17546 16863 23598 
25867 31774 24216 20448 12548 15129 11687 11581

Answer 1

@NickCox has presented an interesting strategy (+1). I might consider it more exploratory in nature however, due to the concern that @whuber points out.

Let me suggest another strategy: You could fit a Gaussian finite mixture model. Note that this makes the very strong assumption that your data are drawn from one or more true normals. As both @whuber and @NickCox point out in the comments, without a substantive interpretation of these data—supported by well-established theory—to support this assumption, this strategy should be considered exploratory as well.

First, let's follow @Glen_b's suggestion and look at your data using twice as many bins:

We still see two modes; if anything, they come through more clearly here. (Note also that the kernel density line should be identical, but appears more spread out due to the larger number of bins.)

Now lets fit a Gaussian finite mixture model. In R, you can use the Mclust package to do this:

library(mclust)
x.gmm = Mclust(x)
summary(x.gmm)
# ----------------------------------------------------
# Gaussian finite mixture model fitted by EM algorithm 
# ----------------------------------------------------
#   
# Mclust V (univariate, unequal variance) model with 2 components:
#   
#   log.likelihood   n df       BIC       ICL
#        -1200.874 120  5 -2425.686 -2442.719
# 
# Clustering table:
#  1  2 
# 68 52 

Two normal components optimizes the BIC. For comparison, we can force a one component fit and perform a likelihood ratio test:

x.gmm.1 = Mclust(x, G=1)
logLik(x.gmm.1)
# 'log Lik.' -1226.241 (df=2)
logLik(x.gmm)-logLik(x.gmm.1)
# 'log Lik.' 25.36657 (df=5)
1-pchisq(25.36657, df=3)  # [1] 1.294187e-05

This suggests it is extremely unlikely you would find data as far from unimodal as yours if they came from a single true normal distribution.

Some people don't feel comfortable using a parametric test here (although if the assumptions hold, I don't know of any problem). One very broadly applicable technique is to use the Parametric Bootstrap Cross-fitting Method (I describe the algorithm here). We can try applying it to these data:

x.gmm$parameters
# $mean
# 12346.98 23322.06 
# $variance$sigmasq
# [1]  4514863 24582180
x.gmm.1$parameters
# $mean
# [1] 17520.91
# $variance$sigmasq
# [1] 43989870

set.seed(7809)
B = 10000;    x2.d = vector(length=B);    x1.d = vector(length=B)
for(i in 1:B){
  x2      = c(rnorm(68, mean=12346.98, sd=sqrt( 4514863)), 
              rnorm(52, mean=23322.06, sd=sqrt(24582180)) )
  x1      = rnorm( 120, mean=17520.91, sd=sqrt(43989870))
  x2.d[i] = Mclust(x2, G=2)$loglik - Mclust(x2, G=1)$loglik
  x1.d[i] = Mclust(x1, G=2)$loglik - Mclust(x1, G=1)$loglik
}
x2.d = sort(x2.d);  x1.d = sort(x1.d)
summary(x1.d)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# -0.29070 -0.02124  0.41460  0.88760  1.36700 14.01000 
summary(x2.d)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  9.006  23.770  27.500  27.760  31.350  53.500 

The summary statistics, and the kernel density plots for the sampling distributions show several interesting features. The log likelihood for the single component model is rarely greater than that of the two component fit, even when the true data generating process has only a single component, and when it is greater, the amount is trivial. The idea of comparing models that differ in their ability to fit data is one of the motivations behind the PBCM. The two sampling distributions barely overlap at all; only .35% of x2.d are less than the maximum x1.d value. If you selected a two component model if the difference in log likelihood were >9.7, you would incorrectly select the one component model .01% and the two component model .02% of the time. These are highly discriminable. If, on the other hand, you chose to use the one component model as a null hypothesis, your observed result is sufficiently small as not to show up in the empirical sampling distribution in 10,000 iterations. We can use the rule of 3 (see here) to place an upper bound on the p-value, namely, we estimate your p-value is less than .0003. That is, this is highly significant.

This raises the question of why these results diverge so much from your dip test. (To answer your explicit question, your dip test provides no evidence that there are two real modes.) I honestly don't know the dip test, so it's hard to say; it may be underpowered. However, I think the likely answer is that this approach assumes your data are generated by true normal[s]. A Shapiro-Wilk test for your data is highly significant ($p < .000001$), and it is also highly significant for the optimal Box-Cox transformation of your data (the inverse square root; $p < .001$). However, data are never really normal (cf., this famous quote), and the underlying components, should they exist, aren't guaranteed to be perfectly normal either. If you find it reasonable that your data could come from a positively skewed distribution, rather than a normal, this level of bimodality may well be within the typical range of variation, which is what I suspect the dip test is saying.

Answer 2

Following up on the ideas in @Nick's answer and comments, you can see how wide the bandwidth needs to be to just flatten out the secondary mode:

Take this kernel density estimate as the proximal null—the distribution closest to the data yet still consistent with the null hypothesis that it's a sample from a unimodal population—and simulate from it. In the simulated samples the secondary mode doesn't often look so distinct, and you needn't widen the bandwidth as much to flatten it out.

Formalizing this approach leads to the test given in Silverman (1981), "Using kernel density estimates to investigate modality", JRSS B, 43, 1. Schwaiger & Holzmann's silvermantest package implements this test, and also the calibration procedure described by Hall & York (2001), "On the calibration of Silverman's test for multimodality", Statistica Sinica, 11, p 515, which adjusts for asymptotic conservatism. Performing the test on your data with a null hypothesis of unimodality results in p-values of 0.08 without calibration and 0.02 with calibration. I'm not familiar enough with the dip test to guess at why it might differ.

R code:

      # kernel density estimate for x using Sheather-Jones 
          # method to estimate b/w:
    density(x, kernel="gaussian", bw="SJ") -> dens.SJ
      # tweak b/w until mode just disappears:
    density(x, kernel="gaussian", bw=3160) -> prox.null
      # fill matrix with simulated samples from the proximal 
      # null:
    x.sim <- matrix(NA, nrow=length(x), ncol=10)
    for (i in 1:10){
      x.sim[ ,i] <- rnorm(length(x), sample(x, size=length(x), 
                          replace=TRUE), prox.null$bw)
    }
      # perform Silverman test without Hall-York calibration:
    require(silvermantest)
    silverman.test(x, k=1, M=10000, adjust=F)
      # perform Silverman test with Hall-York calibration:
    silverman.test(x, k=1, M=10000, adjust=T)

Answer 3

The things to worry about include:

1. The size of the dataset. It is not tiny, not large.

2. The dependence of what you see on histogram origin and bin width. With only one choice evident, you (and we) have no idea of sensitivity.

3. The dependence of what you see on kernel type and width and whatever other choices are made for you in density estimation. With only one choice evident, you (and we) have no idea of sensitivity.

Elsewhere I have suggested tentatively that credibility of modes is supported (but not established) by a substantive interpretation and by the ability to discern the same modality in other datasets of the same size. (Bigger is better too....)

We can't comment on either of those here. One small handle on repeatability is to compare what you get with bootstrap samples of the same size. Here are the results of a token experiment using Stata, but what you see is arbitrarily limited to Stata's defaults, which themselves are documented as plucked out of the air. I got density estimates for the original data and for 24 bootstrap samples from the same.

The indication (no more, no less) is what I think experienced analysts would just guess any way from your graph. The left-hand mode is highly repeatable and the right-hand is distinctly more fragile.

Note that there is an inevitability about this: as there are fewer data nearer the right-hand mode, it won't always reappear in a bootstrap sample. But this is also the key point.

Note that point 3. above remains untouched. But the results are somewhere between unimodal and bimodal.

For those interested, this is the code:

clear 
set scheme s1color 
set seed 2803 

mat data = (10346, 13698, 13894, 19854, 28066, 26620, 27066, 16658, 9221, 13578, 11483, 10390, 11126, 13487, 15851, 16116, 24102, 30892, 25081, 14067, 10433, 15591, 8639, 10345, 10639, 15796, 14507, 21289, 25444, 26149, 23612, 19671, 12447, 13535, 10667, 11255, 8442, 11546, 15958, 21058, 28088, 23827, 30707, 19653, 12791, 13463, 11465, 12326, 12277, 12769, 18341, 19140, 24590, 28277, 22694, 15489, 11070, 11002, 11579, 9834, 9364, 15128, 15147, 18499, 25134, 32116, 24475, 21952, 10272, 15404, 13079, 10633, 10761, 13714, 16073, 23335, 29822, 26800, 31489, 19780, 12238, 15318, 9646, 11786, 10906, 13056, 17599, 22524, 25057, 28809, 27880, 19912, 12319, 18240, 11934, 10290, 11304, 16092, 15911, 24671, 31081, 27716, 25388, 22665, 10603, 14409, 10736, 9651, 12533, 17546, 16863, 23598, 25867, 31774, 24216, 20448, 12548, 15129, 11687, 11581)
set obs `=colsof(data)' 
gen data = data[1,_n] 

gen index = . 

quietly forval j = 1/24 { 
    replace index = ceil(120 * runiform()) 
    gen data`j' = data[index]
    kdensity data`j' , nograph at(data) gen(xx`j' d`j') 
} 

kdensity data, nograph at(data) gen(xx d) 

local xstuff xtitle(data/1000) xla(10000 "10" 20000 "20" 30000 "30") sort 
local ystuff ysc(r(0 .0001)) yla(none) `ystuff'   

local i = 1 
local colour "orange" 
foreach v of var d d? d?? { 
    line `v' data, lc(`colour') `xstuff'  `ystuff' name(g`i', replace) 
    local colour "gs8" 
    local G `G' g`i' 
    local ++i 
} 

graph combine `G' 

Answer 4

LP Nonparametric Mode Identification (name of the algorithm LPMode, the ref of the paper is given below)

MaxEnt Modes [Red color triangles in the plot]: 12783.36 and 24654.28.

L2 Modes [Green color triangles in the plot]: 13054.70 and 24111.61.

Interesting to note the modal shapes, especially the second one which shows considerable skewness (Traditional Gaussian Mixture model likely to fail here).

Mukhopadhyay, S. (2016) Large-Scale Mode Identification and Data-Driven Sciences. https://arxiv.org/abs/1509.06428

Answer 5

I'll add that you can replace the progress bars, which don't illuminate all that much, with an iteration count instead with this:

    x2.d[i] = Mclust(x2, G=2, verbose=FALSE)$loglik - Mclust(x2, 
                     G=1, verbose=FALSE)$loglik
    x1.d[i] = Mclust(x1, G=2, verbose=FALSE)$loglik - Mclust(x1, 
                     G=1, verbose=FALSE)$loglik
    cat(sprintf("\rIteration %s of %s (%.1f%% complete)", i, B, 
              i/B*100))