---
title: Skew-Normal Censored Mixture Model
url: https://www.emergentmind.com/topics/skew-normal-censored-mixture-model
type: topic
---

# Skew-Normal Censored Mixture Model

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" [2009.10826], 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` [2009.10826].

## 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 [2009.10826].

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, $\boldsymbol\Sigma$ is positive definite dispersion, and $\boldsymbol\lambda\in\mathbb R^p$ is the skewness parameter. When $\boldsymbol\lambda=\mathbf 0$, the model reduces to the multivariate normal distribution [2009.10826].

A useful closed form for the cdf is
\[
F_Y(\mathbf y)=P(Y\le \mathbf y)
=
2\, \Phi_{p+1}\!\left( (\mathbf z^\top,0)^\top;\mathbf 0,\mathbf \Omega \right),
\]
with
\[
\mathbf z=\mathbf y-\boldsymbol\mu,
\qquad
\mathbf\Omega=
\begin{pmatrix}
\boldsymbol\Sigma & -\boldsymbol\delta\\
-\boldsymbol\delta^\top & 1
\end{pmatrix},
\qquad
\boldsymbol\delta=
\frac{\boldsymbol\lambda}{\sqrt{1+\boldsymbol\lambda^\top\boldsymbol\lambda}}.
\]

Because conditioning and truncation are central to estimation, the paper introduces the extended skew-normal distribution $ESN_p(\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda,\tau)$, with density
\[
ESN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda,\tau)
=
\xi^{-1}\,
\phi_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma)\,
\Phi_1\!\left(\tau+\boldsymbol\lambda^\top\boldsymbol\Sigma^{-1/2}(\mathbf y-\boldsymbol\mu)\right),
\]
where
\[
\xi=
\Phi_1\!\left(\frac{\tau}{\sqrt{1+\boldsymbol\lambda^\top\boldsymbol\lambda}}\right).
\]
Setting $\tau=0$ 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 $\mathbb A$, the truncated ESN distribution is defined by
\[
f(\mathbf y\mid \mathbf y\in\mathbb A)
=
\frac{ESN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda,\tau)}
{P(Y\in\mathbb A)}
\mathbf 1_{\mathbb A}(\mathbf y),
\]
and the notation $TESN_p(\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda,\tau;[a,b])$ is used for rectangular truncation. A key normalizing constant is
\[
\mathcal L_p(a,b;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda,\tau)
=
\int_a^b ESN_p(\mathbf y;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda,\tau)\,d\mathbf y.
\]
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 $n$ i.i.d. $p$-variate responses
\[
Y_1,\dots,Y_n \sim SN_p(\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda).
\]
For subject $i$, the observed data are $(V_i,C_i)$, where $C_{ik}=0$ means $Y_{ik}$ is observed exactly and $C_{ik}=1$ means $Y_{ik}$ is interval censored,
\[
V_{1ik}\le Y_{ik}\le V_{2ik}.
\]
The response is partitioned as
\[
Y_i=(Y_i^{o\top},Y_i^{c\top})^\top,
\]
with corresponding partitions of $\boldsymbol\mu_i$, $\boldsymbol\Sigma_i$, and $\boldsymbol\lambda_i$.

Using the ESN conditioning result, the paper shows that $Y_i^o$ follows a skew-normal law and that
\[
Y_i^c\mid Y_i^o=y_i^o \sim ESN_{p_i^c}\big(\boldsymbol\mu_i^{co},\boldsymbol\Sigma_i^{cc.o}, \boldsymbol\lambda_i^{c*}, \tau_i^{co}\big),
\]
where
\[
\boldsymbol\mu_i^{co}
=
\boldsymbol\mu_i^c
+
\boldsymbol\Sigma_i^{co}(\boldsymbol\Sigma_i^{oo})^{-1}(y_i^o-\boldsymbol\mu_i^o),
\]
\[
\boldsymbol\Sigma_i^{cc.o}
=
\boldsymbol\Sigma_i^{cc}
-
\boldsymbol\Sigma_i^{co}(\boldsymbol\Sigma_i^{oo})^{-1}\boldsymbol\Sigma_i^{oc},
\]
\[
\tau_i^{co}
=
\tilde{\boldsymbol\lambda}_i^{o\top}(y_i^o-\boldsymbol\mu_i^o),
\qquad
c_i^{oc}
=
(1+\boldsymbol\lambda_i^{c\top}\boldsymbol\Sigma_i^{cc.o}\boldsymbol\lambda_i^c)^{-1/2}.
\]

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:
\[
L_i(\boldsymbol\theta\mid V_i,C_i)
=
L_{p_i^c}(V^c_{1i},V^c_{2i};
\boldsymbol\mu_i^{co},\boldsymbol\Sigma_i^{cc.o},\boldsymbol\lambda_i^{c*},\tau_i^{co})
\;
SN_{p_i^o}(y_i^o;\boldsymbol\mu_i^o,\boldsymbol\Sigma_i^{oo}, c_i^{oc}\boldsymbol\Sigma_i^{oo\,1/2}\tilde{\boldsymbol\lambda}_i^o).
\]
Hence the full log-likelihood is
\[
\ell(\boldsymbol\theta\mid V,C)=\sum_{i=1}^n \log L_i(\boldsymbol\theta\mid V_i,C_i).
\]

Missingness is handled by taking
\[
V_{1ik}=-\infty,\qquad V_{2ik}=+\infty.
\]
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_c^\top,Y_m^\top)^\top$, where $Y_m$ are truly missing values and $Y_c$ are actually censored values. Then $Y_m\mid Y_c,Y_{\text{obs}}$ is non-truncated ESN, while only $Y_c\mid Y_{\text{obs}}$ 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 $T_i$:
\[
Y_i\mid T_i=t_i \sim N_p(\boldsymbol\mu+\boldsymbol\delta t_i,\boldsymbol\Sigma),
\qquad
T_i\sim HN(0,1).
\]
This turns the skew-normal likelihood into a Gaussian complete-data likelihood and enables the EM algorithm [2009.10826].

For the single-component MSNC model, the complete-data parameterization is
\[
\boldsymbol\eta=(\boldsymbol\mu^\top,\boldsymbol\delta^\top,\boldsymbol\Sigma_\Gamma^\top)^\top,
\]
with $\boldsymbol\Sigma_\Gamma=\mathrm{vech}(\boldsymbol\Sigma)$. The complete-data log-likelihood for subject $i$ is
\[
\ell_{ic}(\boldsymbol\eta)
=
-\frac12
\left\{
\log|\boldsymbol\Sigma|
+
(y_i-\boldsymbol\mu-\boldsymbol\delta t_i)^\top
\boldsymbol\Sigma^{-1}
(y_i-\boldsymbol\mu-\boldsymbol\delta t_i)
\right\}
+c.
\]

At iteration $k$, the E-step computes
\[
Q(\boldsymbol\eta\mid \widehat{\boldsymbol\eta}^{(k)})
=
E\big[\ell_c(\boldsymbol\eta)\mid V,C,\widehat{\boldsymbol\eta}^{(k)}\big].
\]
The sufficient conditional expectations are
\[
\widehat y_i^{(k)}=E(Y_i\mid V_i,C_i,\widehat{\boldsymbol\eta}^{(k)}),
\]
\[
\widehat y_i^{2(k)}=E(Y_iY_i^\top\mid V_i,C_i,\widehat{\boldsymbol\eta}^{(k)}),
\]
\[
\widehat t_i^{(k)}=E(T_i\mid V_i,C_i,\widehat{\boldsymbol\eta}^{(k)}),
\qquad
\widehat t_i^{2(k)}=E(T_i^2\mid V_i,C_i,\widehat{\boldsymbol\eta}^{(k)}),
\]
\[
\widehat{ty}_i^{(k)}=E(T_iY_i\mid V_i,C_i,\widehat{\boldsymbol\eta}^{(k)}).
\]
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 $W_{0i}\sim TN_p(\cdot;[v_{1i},v_{2i}])$ and truncated ESN variables $W_i\sim TESN_p(\cdot;[v_{1i},v_{2i}])$. 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 $Y\sim TESN_p(\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda,\tau;[a,b])$ is
\[
E\!\left[
g(Y)\,
\frac{\phi_1(\tau+\boldsymbol\lambda^\top\boldsymbol\Sigma^{-1/2}(Y-\boldsymbol\mu))}
{\Phi_1(\tau+\boldsymbol\lambda^\top\boldsymbol\Sigma^{-1/2}(Y-\boldsymbol\mu))}
\right]
=
\eta\,
\frac{L_p(a,b;-\mathbf b,\mathbf A)}
{\mathcal L_p(a,b;\boldsymbol\mu,\boldsymbol\Sigma,\boldsymbol\lambda,\tau)}
E[g(W)],
\]
with $W\sim TN_p(-\mathbf b,\mathbf A;[a,b])$ and constants $\eta,\mathbf b,\mathbf A$ defined in the paper. These identities are the key to obtaining closed forms for the E-step.

The M-step updates are closed form:
\[
\widehat{\boldsymbol\mu}^{(k+1)}
=
\frac{1}{n}\sum_{i=1}^n
\left(
\widehat y_i^{(k)}-\widehat t_i^{(k)}\widehat{\boldsymbol\delta}^{(k)}
\right),
\]
\[
\widehat{\boldsymbol\delta}^{(k+1)}
=
\left(\sum_{i=1}^n \widehat t_i^{2(k)}\right)^{-1}
\sum_{i=1}^n
\left(
\widehat{ty}_i^{(k)}-\widehat t_i^{(k)}\widehat{\boldsymbol\mu}^{(k+1)}
\right),
\]
\[
\widehat{\boldsymbol\Sigma}^{(k+1)}
=
\frac1n\sum_{i=1}^n
\Big[
\widehat y_i^{2(k)}
+\widehat{\boldsymbol\mu}^{(k)}\widehat{\boldsymbol\mu}^{(k)\top}
+\widehat t_i^{2(k)}\widehat{\boldsymbol\delta}^{(k)}\widehat{\boldsymbol\delta}^{(k)\top}
-\widehat{\boldsymbol\mu}^{(k)}\widehat y_i^{(k)\top}
-\widehat y_i^{(k)}\widehat{\boldsymbol\mu}^{(k)\top}
-\widehat{ty}_i^{(k)}\widehat{\boldsymbol\delta}^{(k)\top}
-\widehat{\boldsymbol\delta}^{(k)}\widehat{ty}_i^{(k)\top}
+\widehat t_i^{(k)}\widehat{\boldsymbol\delta}^{(k)}\widehat{\boldsymbol\mu}^{(k)\top}
+\widehat t_i^{(k)}\widehat{\boldsymbol\mu}^{(k)}\widehat{\boldsymbol\delta}^{(k)\top}
\Big].
\]

Iteration continues until
\[
\left|
\frac{\ell(\widehat{\boldsymbol\eta}^{(k+1)}\mid V,C)}
{\ell(\widehat{\boldsymbol\eta}^{(k)}\mid V,C)}
-1
\right|<\varepsilon,
\]
with a typical choice $\varepsilon=10^{-6}$. 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
\[
Y_i \sim \sum_{j=1}^G \pi_j\, SN_p(\boldsymbol\mu_j,\boldsymbol\Sigma_j,\boldsymbol\lambda_j),
\qquad
\pi_j>0,\ \sum_{j=1}^G\pi_j=1.
\]
Introducing latent class indicators
\[
Z_{ij}=
\begin{cases}
1,& \text{if observation }i\text{ comes from component }j,\\
0,& \text{otherwise},
\end{cases}
\]
the density for the observed censored or missing data becomes
\[
f(V_i\mid C_i,\boldsymbol\theta)
=
\sum_{j=1}^G \pi_j\, f_{ij}(V_i\mid C_i,\boldsymbol\theta),
\]
where $f_{ij}$ is the component-specific MSNC contribution. This is the FM-MSNC model [2009.10826].

The complete data now include $Y_i$, $T_i$, and $Z_{ij}$, with complete-data log-likelihood
\[
\ell_c(\boldsymbol\theta)
=
\sum_{i=1}^n\sum_{j=1}^G Z_{ij}
\left[
\log\pi_j
-\frac12\log|\boldsymbol\Sigma_j|
-\frac12(y_i-\boldsymbol\mu_j-\boldsymbol\delta_j t_i)^\top
\boldsymbol\Sigma_j^{-1}
(y_i-\boldsymbol\mu_j-\boldsymbol\delta_j t_i)
\right]
+c.
\]

The E-step computes posterior component weights
\[
\mathcal Z_{ij}^{(k)}
=
E(Z_{ij}\mid V_i,C_i,\widehat{\boldsymbol\theta}^{(k)})
=
\frac{\widehat\pi_j^{(k)} f_{ij}(V_i\mid C_i,\widehat{\boldsymbol\theta}_j^{(k)})}
{\sum_{\ell=1}^G \widehat\pi_\ell^{(k)} f_{i\ell}(V_i\mid C_i,\widehat{\boldsymbol\theta}_\ell^{(k)})},
\]
together with weighted conditional moments
\[
\mathcal E_{1ij}=E(Z_{ij}Y_i\mid V_i,C_i,\widehat{\boldsymbol\theta}^{(k)}),\quad
\mathcal E_{2ij}=E(Z_{ij}Y_iY_i^\top\mid\cdot),
\]
\[
\mathcal E_{3ij}=E(Z_{ij}T_iY_i\mid\cdot),\quad
\mathcal E_{4ij}=E(Z_{ij}T_i^2\mid\cdot),\quad
\mathcal E_{5ij}=E(Z_{ij}T_i\mid\cdot),
\]
with
\[
\mathcal E_{sij}^{(k)}
=
\mathcal Z_{ij}^{(k)}
\times
E(\cdot \mid V_i,C_i,\widehat{\boldsymbol\theta}^{(k)}, Z_{ij}=1).
\]

The M-step has closed-form updates for $\pi_j$, $\boldsymbol\mu_j$, $\boldsymbol\delta_j$, and $\boldsymbol\Sigma_j$. 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 $\pi_j^{(0)}$, and initialize $\boldsymbol\mu_j,\boldsymbol\Sigma_j,\boldsymbol\lambda_j$ 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
\[
I_e(\widehat{\boldsymbol\theta}\mid y)
=
\sum_{i=1}^n \widehat s_i\,\widehat s_i^\top,
\]
since at the MLE the summed score is zero. Here $\widehat s_i$ is the observed-data score for subject $i$, obtained as the conditional expectation of the complete-data score:
\[
s(y_i\mid \boldsymbol\theta)
=
E\!\left(
\left.
\frac{\partial \ell_i(\boldsymbol\theta\mid y_{ci})}{\partial \boldsymbol\theta}
\right|V_i,C_i,\boldsymbol\theta
\right).
\]
The paper derives explicit score components for mixing weights, component means, covariance parameters via a square-root factor $F_j$, and skewness parameters. Standard errors are obtained from the inverse of $I_e$.

Model selection is based on
- AIC: $-2\ell(\widehat\theta)+2\rho$,
- BIC: $-2\ell(\widehat\theta)+\rho\log n$,
- EDC: $-2\ell(\widehat\theta)+\rho c_n$, with $c_n=0.2\sqrt n$.

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

## 6. Empirical behavior, applications, and implementation

The empirical study in the paper includes simulation experiments and a real-data application [2009.10826]. In simulations with two-component mixtures and left-censoring levels $5\%$ and $30\%$, 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 $n$ reduced bias and variability. Under missing data rates of $5\%$ or $20\%$, 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 $n=500$ with $5\%$ missingness, FM-MSNC gave MAE $=1.9009$ and MARE $=0.8642$, whereas mean imputation gave MAE $=3.2550$ and MARE $=1.6897$.

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 $2$–$4$ 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 $n=300$, FM-MSNC ACCR was about $0.9733$ with no censoring and $0.9545$ at $20\%$ censoring, while FM-MNC was lower, around $0.9248$ to $0.9155$. 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 $n=184$ freshwater streams with $p=5$ 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 $1$–$3$ components were fitted, and all information criteria favored FM-MSNC with $3$ components. The best FM-MSNC model had log-likelihood $-697.68$, AIC $1489.36$, BIC $1640.47$, and EDC $1522.87$. Its estimated mixture proportions were
\[
(\pi_1,\pi_2,\pi_3)=(0.4699,\,0.3440,\,0.1861).
\]
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.

Source: https://www.emergentmind.com/topics/skew-normal-censored-mixture-model