Skew-Normal Censored Mixture Model

• The Skew-Normal Censored Mixture Model is a finite mixture framework that models heterogeneous, asymmetric, and partially observed data using multivariate skew-normal distributions.
• It employs an EM algorithm with latent half-normal augmentation and truncated moment computations to efficiently estimate parameters under censoring and missingness.
• The method is implemented in the R package 'CensMFM' and has shown robust performance in clustering and imputation compared to traditional approaches.

A skew-normal censored mixture model is a finite mixture model for multivariate censored and/or missing data built on the multivariate skew-normal distribution. In the formulation developed in "Finite mixture modeling of censored and missing data using the multivariate skew-normal distribution" (Alencar et al., 2020), the model is designed for data that are simultaneously heterogeneous, asymmetric, and partially observed. It accommodates multimodality via multiple mixture components, skewness via component-specific skewness vectors, and censoring or missingness via interval-censoring likelihood contributions. Missing values are treated as a special case of interval censoring by setting lower and upper bounds to $-\infty$ and $+\infty$, respectively. The resulting framework is fully likelihood-based, uses an analytically simple, yet efficient, EM-type algorithm, and is implemented in the R package CensMFM (Alencar et al., 2020).

1. Statistical setting and motivation

Finite mixture models have been widely used to model and analyze data from heterogeneous populations. In the setting considered here, additional complications arise because data can be missing or subject to upper and/or lower detection limits, and because measures of each population may depart significantly from normality through asymmetric behavior (Alencar et al., 2020).

The proposed skew-normal censored mixture model addresses three features jointly. First, heterogeneity is represented by several latent subpopulations or components. Second, asymmetry is represented by skew-normal component distributions rather than Gaussian components. Third, censoring and missingness are represented through interval-censoring likelihood contributions. A Gaussian mixture can handle heterogeneity but not skewness well; a Student-$t$ mixture can handle outliers or heavy tails but still does not model skewness directly. The finite mixture of multivariate skew-normal components is therefore positioned as a model for skewed heterogeneous populations with censoring or missingness.

The paper also distinguishes between single-component and mixture formulations. For a single population, the model is referred to as MSNC, a model for censored and missing responses under a multivariate skew-normal law. Its finite-mixture extension is called FM-MSNC. This suggests that the model class is intended both for direct density estimation under censoring and for latent-class or clustering analyses when several subpopulations are present.

2. Multivariate skew-normal and extended skew-normal structure

For a $p$-variate random vector $Y$, the multivariate skew-normal density is written as

$$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$

where $\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$ is the $p$-variate normal density, $\Phi_1(\cdot)$ is the standard univariate normal cdf, $\boldsymbol\mu\in\mathbb R^p$ is location, $+\infty$0 is positive definite dispersion, and $+\infty$1 is the skewness parameter. When $+\infty$2, the model reduces to the multivariate normal distribution (Alencar et al., 2020).

A useful closed form for the cdf is

$+\infty$3

with

$+\infty$4

Because conditioning and truncation are central to estimation, the paper introduces the extended skew-normal distribution $+\infty$5, with density

$+\infty$6

where

$+\infty$7

Setting $+\infty$8 gives the usual skew-normal model. The ESN is important because it has closure properties under conditioning, which the skew-normal alone lacks.

For a truncation region $+\infty$9, the truncated ESN distribution is defined by

$t$0

and the notation $t$1 is used for rectangular truncation. A key normalizing constant is

$t$2

The E-step of the estimation procedure relies on moments of these truncated distributions.

3. Censoring, missingness, and the MSNC likelihood

In the single-component formulation, the data consist of $t$3 i.i.d. $t$4-variate responses

$t$5

For subject $t$6, the observed data are $t$7, where $t$8 means $t$9 is observed exactly and $p$0 means $p$1 is interval censored,

$p$2

The response is partitioned as

$p$3

with corresponding partitions of $p$4, $p$5, and $p$6.

Using the ESN conditioning result, the paper shows that $p$7 follows a skew-normal law and that

$p$8

where

$p$9

$Y$0

$Y$1

The individual contribution to the likelihood is the product of a truncated ESN probability for the censored part conditional on the observed part and the skew-normal density for the observed part: $Y$2 Hence the full log-likelihood is

$Y$3

Missingness is handled by taking

$Y$4

The paper points out, however, that if missing entries are treated as censoring, the dimension of the truncated integrals can become unnecessarily large. To reduce computational burden, it partitions the censored vector as $Y$5, where $Y$6 are truly missing values and $Y$7 are actually censored values. Then $Y$8 is non-truncated ESN, while only $Y$9 requires truncated moments. A plausible implication is that the computational bottleneck is tied primarily to genuinely censored coordinates rather than to all incompletely observed coordinates.

4. Latent-variable construction and EM estimation

The skew-normal model is represented through a latent half-normal variable $$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$0: $$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$1 This turns the skew-normal likelihood into a Gaussian complete-data likelihood and enables the EM algorithm (Alencar et al., 2020).

For the single-component MSNC model, the complete-data parameterization is

$$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$2

with $$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$3. The complete-data log-likelihood for subject $$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$4 is

$$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$5

At iteration $$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$6, the E-step computes

$$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$7

The sufficient conditional expectations are

$$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$8

$$SN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda) = 2\,\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\, \Phi_1\!\left( \boldsymbol\lambda^\top \boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu) \right),$$9

$\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$0

$\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$1

These are obtained using propositions on truncated ESN moments.

The E-step distinguishes three cases: no censored components, all components censored, and mixed observed and censored components. In computational terms, the major point is that the E-step reduces to evaluating the first two moments of truncated multivariate normal variables $\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$2 and truncated ESN variables $\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$3. The paper explicitly uses results from Galarza et al. for these moments, implemented in the R package MomTrunc.

One of the identities used for a truncated ESN variable $\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$4 is

$\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$5

with $\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$6 and constants $\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$7 defined in the paper. These identities are the key to obtaining closed forms for the E-step.

The M-step updates are closed form: $\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$8

$\phi_p(\cdot;\boldsymbol\mu,\boldsymbol\Sigma)$9

$p$0

Iteration continues until

$p$1

with a typical choice $p$2. The working parameters are then transformed back to the original skew-normal parameters through the recovery formulas given in the paper.

5. Finite-mixture extension and inferential framework

The finite-mixture model is

$p$3

Introducing latent class indicators

$p$4

the density for the observed censored or missing data becomes

$p$5

where $p$6 is the component-specific MSNC contribution. This is the FM-MSNC model (Alencar et al., 2020).

The complete data now include $p$7, $p$8, and $p$9, with complete-data log-likelihood

$\Phi_1(\cdot)$0

The E-step computes posterior component weights

$\Phi_1(\cdot)$1

together with weighted conditional moments

$\Phi_1(\cdot)$2

$\Phi_1(\cdot)$3

with

$\Phi_1(\cdot)$4

The M-step has closed-form updates for $\Phi_1(\cdot)$5, $\Phi_1(\cdot)$6, $\Phi_1(\cdot)$7, and $\Phi_1(\cdot)$8. The paper also warns that mixture likelihoods can be unbounded, so the EM solution is treated as the ML estimate in the usual practical sense, with suitable initialization. The recommended initialization is: cluster the data using K-means, use cluster proportions as $\Phi_1(\cdot)$9, and initialize $\boldsymbol\mu\in\mathbb R^p$0 using the mixsmsn R package.

For asymptotic covariance, the paper uses an information-based empirical matrix. Following Meilijson/Basford and Louis’ missing-data score identity, the empirical observed information is approximated by

$\boldsymbol\mu\in\mathbb R^p$1

since at the MLE the summed score is zero. Here $\boldsymbol\mu\in\mathbb R^p$2 is the observed-data score for subject $\boldsymbol\mu\in\mathbb R^p$3, obtained as the conditional expectation of the complete-data score: $\boldsymbol\mu\in\mathbb R^p$4 The paper derives explicit score components for mixing weights, component means, covariance parameters via a square-root factor $\boldsymbol\mu\in\mathbb R^p$5, and skewness parameters. Standard errors are obtained from the inverse of $\boldsymbol\mu\in\mathbb R^p$6.

Model selection is based on

• AIC: $\boldsymbol\mu\in\mathbb R^p$7,
• BIC: $\boldsymbol\mu\in\mathbb R^p$8,
• EDC: $\boldsymbol\mu\in\mathbb R^p$9, with $+\infty$00.

These criteria are used both for selecting the number of mixture components and for comparing FM-MSNC against normal and $+\infty$01-mixture competitors.

6. Empirical behavior, applications, and implementation

The empirical study in the paper includes simulation experiments and a real-data application (Alencar et al., 2020). In simulations with two-component mixtures and left-censoring levels $+\infty$02 and $+\infty$03, the EM estimates were close to the truth, and the approximate standard errors from the information matrix matched Monte Carlo SDs reasonably well. Increasing censoring increased bias somewhat, whereas larger $+\infty$04 reduced bias and variability. Under missing data rates of $+\infty$05 or $+\infty$06, estimates were stable, the standard-error approximation worked well, and skewness parameters were more affected by missingness.

For prediction of missing values, FM-MSNC imputation was compared against mean imputation using MAE and MARE. FM-MSNC was substantially better in all settings. At $+\infty$07 with $+\infty$08 missingness, FM-MSNC gave MAE $+\infty$09 and MARE $+\infty$10, whereas mean imputation gave MAE $+\infty$11 and MARE $+\infty$12.

In selecting the number of components, simulated data generated from a two-component FM-MSNC model led AIC, BIC, and EDC to almost always select the correct FM-MSNC(2) model over normal mixtures with $+\infty$13–$+\infty$14 components. For clustering, correct classification rates were high for FM-MSNC and better than FM-MNC. ACCR decreased as censoring increased but improved with larger sample size. At $+\infty$15, FM-MSNC ACCR was about $+\infty$16 with no censoring and $+\infty$17 at $+\infty$18 censoring, while FM-MNC was lower, around $+\infty$19 to $+\infty$20. Bias and MSE decreased with increasing sample size, supporting consistency and asymptotic validity of the EM estimators.

The real-data application analyzes a dataset of $+\infty$21 freshwater streams with $+\infty$22 trace metals: Cu, Pb, Zn, Ca, and Mg. The data have substantial left censoring, especially for Pb and Zn, and were standardized. FM-MSNC, FM-MNC, and FM-MtC models with $+\infty$23–$+\infty$24 components were fitted, and all information criteria favored FM-MSNC with $+\infty$25 components. The best FM-MSNC model had log-likelihood $+\infty$26, AIC $+\infty$27, BIC $+\infty$28, and EDC $+\infty$29. Its estimated mixture proportions were

$+\infty$30

The fitted clusters revealed strong heterogeneity in the metal concentration profiles, and the paper reports scatterplots and histograms of the predicted clusters. The authors note that clustering is visually difficult in raw data, but the fitted FM-MSNC model separates the subgroups well.

The proposed methods are implemented in the R package CensMFM, and truncated-moment computations rely on MomTrunc. In practical terms, the model is positioned as a tool for censored data, interval-censored data, missing data, and mixture modeling with skewness. The principal methodological contribution is a fully likelihood-based finite mixture of multivariate skew-normal components for censored and missing data with interval-censoring formulation, latent half-normal augmentation, closed-form M-step updates, truncated skew-normal moment calculations in the E-step, information-based standard errors, and software support.