--- title: "Introduction to DEmixR" author: "Farrokh Habibzadeh" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction to DEmixR} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE, fig.width = 6, fig.height = 4 ) ``` # 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. 1. **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. 2. **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 ```{r, eval = FALSE} 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`. ```{r} library(DEmixR) set.seed(123) x <- c(rnorm(150, mean = 0, sd = 1), rnorm(100, mean = 4, sd = 1)) mix2_bounds(x, "normal") ``` 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. ```{r} 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"*): ```{r} 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. ```{r} sel <- select_best_mixture(x, n_runs = 1, NP = 25, itermax = 300, quiet = 0) sel # uses print.demixr_select sel$best$family ``` # Fitting a normal mixture ```{r} fit <- fit_norm2(x, n_runs = 1, NP = 25, itermax = 300, quiet = 0) fit # print.demixr_fit: compact overview summary(fit) # summary.demixr_fit: per-component table 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. ```{r} plot(fit, which = "pit") plot(fit, which = "qq") ``` The pieces are also available directly if you need them: ```{r} fit$par # named vector: p, m1, s1, m2, s2 fit$logLik # exact log-likelihood (log-sum-exp) c(AIC = fit$AIC, BIC = fit$BIC) fit$source # "DEoptim" = converged L-BFGS-B refinement of the DE solution fit$de_logLik # log-likelihood reached by each DE run fit$diagnostics$comp_n # expected observations per component ``` # Fitting a lognormal mixture ```{r} 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 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. ```{r} set.seed(123) boot <- bootstrap_mix2(fit_ln, B = 40, parametric = TRUE, quiet = 0) boot # print.demixr_boot ``` # 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)`. ```{r} ev <- evaluate_init(par_init = c(0.5, 0, 0.5, 1.6, 0.4), x = y, family = "lognormal") ev$success ev$logLik ``` # 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). ```r 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. ```r 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 ```r # 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. # Limitations - Exactly two components; normal and lognormal families only. - The estimate is the best solution inside the bounds, not an unrestricted maximum-likelihood estimate; sensitivity to `scale_range` and `p_range` should be checked. - Differential evolution is stochastic and has no optimality guarantee; it is slower than a single EM fit. - Bootstrap confidence intervals are percentile intervals; with `refit = "local"` they do not reflect multimodality. - Missing values are not accepted (`x` must be finite). # References 1. Hathaway, R. J. (1985). A constrained formulation of maximum-likelihood estimation for normal mixture distributions. *The Annals of Statistics*, 13(2), 795--800. For comparisons of DEmixR with EM, multi-start EM (Biernacki, Celeux & Govaert, 2003, *Computational Statistics & Data Analysis*, 41, 561--575) and the benchmark of Gaussian-mixture software in Chassagnol et al. (2023, *The R Journal*, 15(2), 56--76) are the relevant references. 2. 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. 3. R Core Team (2024). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria.