Introduction

DEmixR fits and evaluates two-component mixture models (normal and lognormal). Rather than the usual Expectation-Maximization (EM) algorithm, it searches the parameter space with a global optimizer — Differential Evolution, via the DEoptim package — and then refines the best candidate with a local quasi-Newton ("L-BFGS-B") step. Repeating the global search from several random populations (n_runs) makes the fit less dependent on a single set of starting values.

Two things are worth knowing from the outset.

The search is bounded. The likelihood of a normal or lognormal mixture with component-specific variances is unbounded (a component can collapse onto one observation), so there is no finite global maximum of the unrestricted likelihood. DEmixR therefore returns the best solution found inside a box of data-derived bounds (see Model and parameter bounds). It is a regularized estimate, not an unrestricted maximum-likelihood estimate.
Global search is not a guarantee, and it is not the only remedy for local optima. A multi-start EM algorithm run with a comparable computational budget is a legitimate alternative, and the general mixture packages mixtools, mclust, mixR and flexmix cover more components and distributions. DEmixR is deliberately specialized.

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

Installation

install.packages("DEmixR")

Overview of the exported functions

Function	Purpose
`prelim_plots()`	Exploratory plots (histogram, normal Q-Q and P-P, log-scale versions)
`mix2_bounds()`	Show (and modify) the parameter bounds used by the fit
`select_best_mixture()`	Compare lognormal vs normal mixtures by BIC (or AIC)
`fit_norm2()`	Fit a two-component normal mixture
`fit_lognorm2()`	Fit a two-component lognormal mixture
`bootstrap_mix2()`	Bootstrap the fitted parameters
`evaluate_init()`	Refine a specific set of starting values

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 (plot(fit, which = "density" | "pit" | "qq")).

Model and parameter bounds

For \(x\) (normal family) with sample standard deviation \(s_x\) the search is carried out in the box returned by mix2_bounds():

mixing proportion \(p \in [0.001, 0.999]\);
means \(m_1 \in [\min x - 2s_x, \max x]\) and \(m_2 \in [\min x, \max x + 2s_x]\) (the means are not restricted to opposite sides of the median; labels are assigned afterwards so that \(m_1 \le m_2\));
standard deviations \(s_1, s_2 \in [0.1 s_x, 10 s_x]\).

For the lognormal family the same rule is applied to \(\log x\), and \(m_j\), \(s_j\) are the meanlog and sdlog of component \(j\). The factors 0.1 and 10 are pragmatic regularization choices; they can be changed with scale_range and p_range, or the whole box can be supplied through lower and upper.

library(DEmixR)
set.seed(123)
x <- c(rnorm(150, mean = 0, sd = 1),
       rnorm(100, mean = 4, sd = 1))
mix2_bounds(x, "normal")
#> $lower
#>          p         m1         s1         m2         s2 
#>  0.0010000 -6.6996725  0.2195252 -2.3091689  0.2195252 
#> 
#> $upper
#>        p       m1       s1       m2       s2 
#>  0.99900  7.24104 21.95252 11.63154 21.95252

Optionally, a relative constraint on the scale ratio in the sense of Hathaway (1985), scale_ratio, requires \(\min(s_1, s_2)/\max(s_1, s_2)\) to be at least the given value (default 0 = off). It removes very narrow components next to wide ones, but it also excludes any true solution whose scale ratio is smaller, so it should be chosen with the subject matter in mind.

A fitted solution that lies on one of these bounds was determined by the bound, not by the data; the fitting functions warn when that happens, and when a component’s expected size (the sum of its posterior probabilities) is below min_comp_n = 5 observations. With small samples and a minor component, the best solution inside the box can be a spurious component that sits on a few extreme observations.

Diagnostic plots

Start by looking at the data. These plots compare the data with a single normal (or, on the log scale, lognormal) distribution. A mixture is not itself normal, so curvature in the Q-Q or P-P plot is expected; the plots are exploratory and do not validate the assumed component family.

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 (or AIC, with criterion = "AIC"). It uses the same search settings as the fitting functions. When the data contain non-positive values the lognormal family is impossible, so the choice is then dictated by the support of the data rather than by the criterion. Its print() method reports the criterion 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
#> (lognormal family not considered: x contains non-positive values)
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
#> Likelihood: exact; best solution found within the parameter bounds
#> 
#> 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
#> Expected observations per component (sum of posterior probabilities): 142.6 / 107.4
plot(fit)           # plot.demixr_fit: fitted mixture over a histogram

Two diagnostics that compare the data with the fitted mixture (rather than with a single normal distribution) are available through which: the histogram of the probability integral transform (approximately uniform when the fit describes the data) and a Q-Q plot against the fitted mixture quantiles.

plot(fit, which = "pit")

plot(fit, which = "qq")

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          # exact log-likelihood (log-sum-exp)
#> [1] -500.1969
c(AIC = fit$AIC, BIC = fit$BIC)
#>      AIC      BIC 
#> 1010.394 1028.001
fit$source          # "DEoptim" = converged L-BFGS-B refinement of the DE solution
#> [1] "DEoptim"
fit$de_logLik       # log-likelihood reached by each DE run
#> [1] -500.1969
fit$diagnostics$comp_n   # expected observations per component
#>    comp1    comp2 
#> 142.6142 107.3858

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
#> Likelihood: exact; best solution found within the parameter bounds
#> 
#> 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 the bootstrap median (or mean, center = "mean") and percentile confidence intervals. (B is kept small here for speed.)

By default (refit = "local") every replicate is refitted by L-BFGS-B started at the original estimate, which is fast but does not repeat the global search, so it does not propagate uncertainty due to multimodality. With refit = "de" every replicate is refitted with the same procedure as the original fit (much slower). Percentile intervals of mixture parameters can perform poorly with weakly separated components or small mixing proportions.

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 converged)
#> Resampling: parametric; refit: local; 95% percentile intervals
#> 
#>             p      m1     s1     m2     s2
#> median 0.5439 -0.1019 0.4117 1.5110 0.4468
#> 2.5%   0.4926 -0.1493 0.3451 1.4171 0.3591
#> 97.5%  0.6117 -0.0380 0.4935 1.6345 0.5397
#> 
#> Note: replicates were refitted locally from the original estimate; the global
#> search was not repeated (see ?bootstrap_mix2, argument 'refit').

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.

Tuning the fit

No seed argument is provided (the differential-evolution search is stochastic); call set.seed() immediately before the fit for a reproducible result (guaranteed for parallelType = 0, the default).

set.seed(1)
fit_ln <- fit_lognorm2(
  y,
  NP = 150,           # DEoptim population size
  n_runs = 20,        # independent runs (default 10)
  itermax = 2000,     # maximum iterations per run
  parallelType = 0,   # 0 = serial; see below
  quiet = 2,          # verbosity (see below)
  par_init = NULL,    # optional starting values for an extra local search
  scale_range = c(0.1, 10)   # bounds of the component sds (multiples of sd(x))
)

Key options

NP: larger populations explore more thoroughly but run more slowly.
n_runs: more runs improve robustness at a roughly proportional time cost.
parallelType: passed to DEoptim::DEoptim.control(). 0 = serial, 1 = the parallel package (needs the parallelly package, which DEoptim only suggests), 2 = foreach (needs a registered backend).
quiet: 0 = silent, 1 = per-run messages plus the DEoptim trace, 2 = a text progress bar.

Bootstrapping

No seed argument is provided here either; call set.seed() before the call for reproducible resampling.

set.seed(1)
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
  refit = "de",       # repeat a (short) global search for every replicate
  de_control = list(NP = 50, n_runs = 1, itermax = 2000),
  parallelType = 4,   # integer > 1: number of forked worker processes (Unix)
  quiet = 2
)

boot$central          # median (default) or mean of the bootstrap distribution
boot$ci               # percentile confidence intervals
boot$n_failed         # replicates without a converged refinement

Checking starting values

# 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 (code 0)?
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 global search reduces dependence on a single set of starting values, but it does not remove multimodality, it is carried out inside data-derived bounds, and it is not necessarily better than a multi-start EM algorithm run for a comparable time. The two components of an unlabeled mixture are statistical components; interpreting them as, for example, healthy and diseased subpopulations requires external validation.

Introduction to DEmixR

Farrokh Habibzadeh

2026-09-25