Introduction

DEmixR fits and evaluates two-component mixture models (normal and lognormal). Rather than the usual Expectation-Maximization (EM) algorithm, it locates the maximum-likelihood estimates with a global search — Differential Evolution, via the DEoptim package — followed by a local quasi-Newton ("L-BFGS-B") refinement. Global search makes the fit less sensitive to starting values and local optima, which can be helpful when the two components overlap strongly.

The code below uses deliberately small data sets and small optimizer settings so that the vignette builds quickly. For real analyses, increase n_runs, NP, and itermax (see Advanced usage).

Installation

install.packages("DEmixR")

Overview of the exported functions

The objects returned by fit_norm2() / fit_lognorm2() (class demixr_fit), bootstrap_mix2() (class demixr_boot) and select_best_mixture() (class demixr_select) all carry print() methods, and a demixr_fit additionally has summary() and plot() methods.

library(DEmixR)

Diagnostic plots

Start by looking at the data.

set.seed(123)
x <- c(rnorm(150, mean = 0, sd = 1),
       rnorm(100, mean = 4, sd = 1))

prelim_plots(x, which = c("hist", "qq"))

Passing col_density = NA suppresses the kernel-density overlay (and the histogram title changes from “Histogram with density” to “Histogram”):

prelim_plots(x, which = "hist", col_density = NA)

Model selection

select_best_mixture() fits both families and keeps the one with the lower BIC. Its print() method reports the BIC of each family and the winner.

sel <- select_best_mixture(x, n_runs = 1, NP = 25, itermax = 300, quiet = 0)
sel                 # uses print.demixr_select
#> DEmixR model selection
#> 
#>   normal 
#> 1028.001 
#> 
#> Preferred family (lowest BIC): normal
sel$best$family
#> [1] "normal"

Fitting a normal mixture

fit <- fit_norm2(x, n_runs = 1, NP = 25, itermax = 300, quiet = 0)

fit                 # print.demixr_fit: compact overview
#> Two-component normal mixture (DEmixR)
#> Convergence: successful
#> 
#> Parameter estimates:
#>       p      m1      s1      m2      s2 
#>  0.5705 -0.1173  0.8715  3.8607  1.0668 
#> 
#> logLik = -500.1969   AIC = 1010.3939   BIC = 1028.0012
summary(fit)        # summary.demixr_fit: per-component table
#> Two-component normal mixture fitted to n = 250 observations
#> 
#>  component weight location  scale
#>          1 0.5705  -0.1173 0.8715
#>          2 0.4295   3.8607 1.0668
#> 
#> logLik = -500.1969, AIC = 1010.3939, BIC = 1028.0012
plot(fit)           # plot.demixr_fit: fitted mixture over a histogram

The pieces are also available directly if you need them:

fit$par             # named vector: p, m1, s1, m2, s2
#>          p         m1         s1         m2         s2 
#>  0.5704521 -0.1173263  0.8714674  3.8607499  1.0667846
fit$logLik
#> [1] -500.1969
c(AIC = fit$AIC, BIC = fit$BIC)
#>      AIC      BIC 
#> 1010.394 1028.001

Fitting a lognormal mixture

set.seed(123)
y <- c(rlnorm(150, meanlog = 0, sdlog = 0.5),
       rlnorm(100, meanlog = 1.6, sdlog = 0.4))

fit_ln <- fit_lognorm2(y, n_runs = 1, NP = 25, itermax = 300, quiet = 0)
fit_ln
#> Two-component lognormal mixture (DEmixR)
#> Convergence: successful
#> 
#> Parameter estimates:
#>       p      m1      s1      m2      s2 
#>  0.5491 -0.0879  0.4157  1.5156  0.4387 
#> 
#> logLik = -453.7858   AIC = 917.5716   BIC = 935.1789
plot(fit_ln)

Bootstrap confidence intervals

bootstrap_mix2() returns central estimates and percentile confidence intervals; its print() method shows both. (B is kept small here for speed.)

set.seed(123)
boot <- bootstrap_mix2(fit_ln, B = 40, parametric = TRUE, quiet = 0)
boot                # print.demixr_boot
#> Bootstrap for a two-component lognormal mixture (B = 40, 40 successful)
#> 
#>              p      m1     s1     m2     s2
#> central 0.5452 -0.0789 0.4113 1.5044 0.4412
#> 2.5%    0.4599 -0.1929 0.3571 1.4103 0.3585
#> 97.5%   0.6431  0.0201 0.5015 1.6134 0.4936

Evaluating a set of starting values

evaluate_init() runs a single local optimization from starting values you supply — useful for checking whether a particular initialization converges. The parameter vector is c(p, m1, s1, m2, s2).

ev <- evaluate_init(par_init = c(0.5, 0, 0.5, 1.6, 0.4), x = y,
                    family = "lognormal")
ev$success
#> [1] TRUE
ev$logLik
#> [1] -453.7858

Advanced usage

The chunks in this section are shown for reference and are not executed.

fit_ln <- fit_lognorm2(
  y,
  NP = 150,           # DEoptim population size
  n_runs = 20,        # independent runs to reduce the chance of a local optimum
  itermax = 2000,     # maximum iterations per run
  parallelType = 1,   # parallelize the DEoptim runs
  quiet = 2,          # verbosity (see below)
  par_init = NULL      # optional starting values for an extra local search
)

boot <- bootstrap_mix2(
  fit_ln,
  B = 500,            # number of bootstrap replicates
  parametric = TRUE,  # parametric (TRUE) or nonparametric (FALSE) resampling
  ci_level = 0.90,    # confidence level
  parallelType = 1,
  quiet = 2
)

boot$central          # means (normal) or medians (lognormal)
boot$ci               # percentile confidence intervals

# evaluate_init() takes a single starting vector c(p, m1, s1, m2, s2),
# NOT a number of random starts.
ev <- evaluate_init(
  par_init = c(0.5, 0, 1, 4, 1),
  x = x,
  family = "normal"
)

ev$success            # did L-BFGS-B converge?
ev$par                # refined parameters

Summary

DEmixR provides a compact workflow for two-component mixtures: explore the data with prelim_plots(), pick a family with select_best_mixture(), fit it with fit_norm2() / fit_lognorm2(), inspect the result with the print()/summary()/plot() methods, and quantify uncertainty with bootstrap_mix2(). The use of global optimization is intended to make the fit robust to difficult likelihood surfaces.

References

1. Mullen, K. M., Ardia, D., Gil, D. L., Windover, D., & Cline, J. (2011). DEoptim: An R Package for Global Optimization by Differential Evolution. Journal of Statistical Software, 40(6), 1–26. https://doi.org/10.18637/jss.v040.i06
2. R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.