Publications · Conference

A Unified Framework for Variable Selection in Model-Based Clustering with Missing Not at Random

Binh H. Ho⋆, Long Nguyen Chi⋆, TrungTin Nguyen⋆†, Binh T. Nguyen, Van Ha Hoang, Christopher Drovandi

⋆ Co-first author, † Corresponding author.

NeurIPS 2025 · Poster NeurIPS 2025. Advances in Neural Information Processing Systems 38 (2025).

Abstract

Model-based clustering integrated with variable selection is a powerful tool for uncovering latent structures within complex data. However, its effectiveness is often hindered by challenges such as identifying relevant variables that define heterogeneous subgroups and handling data that are missing not at random, a prevalent issue in fields like transcriptomics. While several notable methods have been proposed to address these problems, they typically tackle each issue in isolation, thereby limiting their flexibility and adaptability. This paper introduces a unified framework designed to address these challenges simultaneously. Our approach incorporates a data-driven penalty matrix into penalized clustering to enable more flexible variable selection, along with a mechanism that explicitly models the relationship between missingness and latent class membership. We demonstrate that, under certain regularity conditions, the proposed framework achieves both asymptotic consistency and selection consistency, even in the presence of missing data. This unified strategy significantly enhances the capability and efficiency of model-based clustering, advancing methodologies for identifying informative variables that define homogeneous subgroups in the presence of complex missing data patterns. The performance of the framework, including its computational efficiency, is evaluated through simulations and demonstrated using both synthetic and real-world transcriptomic datasets. ††⋆Co-first author, †Corresponding author.

1 Introduction

Model-based Clustering. Model-based clustering formulates clustering as a probabilistic inference task, assuming data are generated from a finite mixture model, with each component corresponding to a cluster. This enables likelihood-based estimation and principled model selection. Gaussian mixture models (GMMs) [42, 43] are a classical instance, which can be estimated via the expectation-maximization (EM) algorithm [16], providing soft assignments and flexibility in modeling non-spherical and overlapping clusters. Bayesian mixture models extend this framework by treating parameters, and even the number of components, as random variables, thereby enabling the quantification of uncertainty. Techniques such as Markov chain Monte Carlo (MCMC) and variational inference are employed to sample the posterior distribution, although challenges such as label switching persist [64]. Bayesian nonparametric approaches, including reversible-jump MCMC and birth-death processes, impose prior constraints on model complexity and infer the number of clusters directly. Despite offering interpretability and robustness, especially in small-sample settings, Bayesian methods suffer from high computational costs, label switching, and convergence issues [29]. Therefore, this paper focuses on the maximum likelihood estimation (MLE) approach for dealing with variable selection and missing data in model-based clustering.

Related Works on Variable Selection for Model-based Clustering. Traditional variable selection methods like best-subset and stepwise regression rely on information criteria such as Akaike information criterion (AIC) [1] and Bayesian information criterion (BIC) [58], but become computationally infeasible as dimensionality increases. Penalized likelihood approaches, notably least absolute shrinkage and selection operator (LASSO) [70], address scalability by inducing sparsity; other earlier contributions include ridge regression [28] and the nonnegative garrote [9]. Efficient algorithms like least-angle regression [18] and coordinate descent enable application in high dimensions. To reduce dimensionality before modeling, screening methods such as sure independence screening [19] and its conditional extension [5] filter variables based on marginal or conditional correlations. In clustering, Law et al. [32] propose simultaneous variable selection and clustering using feature saliency. Andrews and McNicholas [3] introduce variable selection for clustering and classification, combining filter and wrapper strategies to discard noisy variables. Bayesian methods such as those by Tadesse et al. [67] and Kim et al. [30] extend variable relevance indicators to mixture models, though MCMC scalability remains a challenge. Raftery and Dean [54] develop a BIC-based stepwise method for identifying clustering-relevant variables, later optimized by Scrucca [59]. Finally, Maugis et al. [38] enrich this framework by incorporating redundancy modeling, enhancing selection accuracy in correlated settings at the cost of increased computational complexity. To overcome this computational complexity, Celeux et al. [13] propose a two-step strategy for selecting variables in mixture models. First, variables are ranked by optimizing a penalized likelihood function that shrinks both component means and precisions, following the approach of Zhou, Pan, and Shen [75]. This method is scalable to moderate dimensions and yields an ordered list where informative variables are expected to appear early. Second, roles are assigned in a single linear pass through the ranked list, replacing combinatorial search with a faster and interpretable procedure that maintains competitive clustering accuracy. However, there is no theoretical guarantee that this ranking recovers the true signal-redundant-uninformative partition.

Handling Missing Data. Missing data are commonly classified into three categories: missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR). The last category, MNAR, is the most challenging, as the probability of missingness depends on unobserved values. Under the MAR assumption, the method of multiple imputation provides a principled framework by replacing each missing value with several plausible alternatives, then combining analyses across imputed datasets to account for uncertainty [57]. The multivariate imputation by chained equations algorithm [71] is a flexible implementation that fits a sequence of conditional models, accommodating mixed data types. In high-dimensional settings, estimating full joint or conditional models becomes unstable. To address this, regularized regression techniques, such as the LASSO and its Bayesian counterpart, have been integrated into imputation models to improve predictive performance by selecting and shrinking predictors [74]. Nonparametric and machine learning-based imputation methods have also been developed. The MissForest algorithm [63] uses random forests to iteratively impute missing values, capturing nonlinear relations. Deep generative models, such as the missing data importance-weighted autoencoder [36] and its federated variant [4], rely on latent variable representations to generate imputations, assuming data lie near a low-dimensional manifold. For MNAR data, model-based approaches face difficulties due to unidentifiability without external constraints. Two common strategies include selection models, which jointly specify the data and missingness mechanism, and pattern mixture models, which condition on observed missingness patterns. Both frameworks require strong assumptions or instruments to yield valid inferences. To address this challenge, Sportisse et al. [62] studied the identifiability of selection models under the MNAR assumption, particularly when missingness depends on unobserved values. They showed that identifiability can be achieved under structural assumptions, such as requiring each variable’s missingness to be explained either by its value or by the latent class, but not both simultaneously. Their proposed joint modeling framework, termed MNARz, ensures that cluster-specific missingness reflects meaningful distributional differences. A key insight is that under MNARz, the missing data problem can be reformulated as MAR on an augmented data matrix, enabling tractable inference.

Main Contribution. Inspired by the work of Maugis et al. [39], who addressed variable selection under MAR without imputation, we generalize their framework to latent class MNARz. Under the MNARz mechanism, we extend the LASSO-based variable selection framework for model-based clustering proposed by [13], establishing identifiability and consistency guarantees in Section 4 for recovering the true signal-redundant-uninformative partition. Our method achieves substantial empirical improvements in Section 5 over existing approaches. This yields a unified, high-dimensional model-based clustering method that jointly addresses variable selection and MNAR inference in a principled manner.

Paper Organization. The remainder of the paper is structured as follows. Section 2 reviews variable selection in model-based clustering with MNAR data. Section 3 and Section 4 present our proposed framework and its theoretical guarantees, respectively. Simulation and real data results are reported in Section 5. We conclude with a summary, limitations, and future directions in Section 6. Proofs and additional details are provided in the supplementary material.

Notation. Throughout the paper, we use the shorthand [N] to denote the index set {1,2,…,N} for any positive integer N∈ℕ. We denote ℝ as the set of real numbers. The notation |𝕊| represents the cardinality of any set 𝕊. For any vector 𝒗∈ℝD, ‖𝒗‖p denotes its p-norm. The operator ⊙ indicates the element-wise (Hadamard) product between two matrices. We use n∈[N] to index observations and d∈[D] to index variables. The vector 𝒚n∈ℝD denotes the full data vector for observation n, which may be further specified in more granular forms. The binary vector 𝒄n=(cn​1,…,cn​D)∈{0,1}D represents the missingness mask, where cn​d=1 if yn​d is missing. The latent variable zn∈[K] denotes the mixture component assignment for observation n, with indicator variable zn​k=𝟏{zn=k}. We use ℓ​(⋅) to denote the log-likelihood.

2 Background

We begin by introducing some preliminaries on model-based clustering with variable selection, penalized clustering, and MNARz formulation.

Model-based Clustering and Parsimonious Mixture Form. Model-based clustering adopts a probabilistic framework by assuming that the data matrix 𝐘=(𝒚1,…,𝒚N)⊤∈ℝN×D, with each 𝒚n∈ℝD, is independently drawn from a location-scale mixture distribution, a class known for its universal approximation capabilities and favorable convergence properties [25, 55, 14, 53, 52, 46, 47, 49, 60, 26, 27]. The density of an observation 𝒚n under a K-component GMM with parameters 𝜶={𝝅,𝝁,𝚺} is given by: fGMM​(𝒚n∣K,m,𝜶)=∑k=1Kπk​𝒩​(𝒚n∣𝝁k,𝚺k) where πk>0 and ∑k=1Kπk=1. The term 𝒩​(𝒚n∣𝝁k,𝚺k) denotes the multivariate normal density with mean vector 𝝁k and covariance matrix 𝚺k, which is compactly denoted by 𝚯k=(𝝁k,𝚺k). The covariance matrix 𝚺k encodes the structure determined by model form m, allowing for parsimonious modeling via spectral decompositions that control cluster volume, shape, and orientation [12]. These constraints are particularly valuable in high-dimensional settings, where D≫N and standard GMMs require estimating 𝒪​(K2​D) parameters, leading to overparameterization. Model selection criteria such as the BIC [58, 39, 22, 21, 48], slope heuristics [6, 51, 50], integrated classification likelihood (ICL) [7, 24], extended BIC (eBIC) [23, 45], dendrogram selection criterion [17, 69], and Sin-White information criterion (SWIC) [61, 72] can be employed to select both the number of components K and the covariance structure m.

Variable Selection as a Model Selection Problem in Model-based Clustering. In high-dimensional settings, many variables may be irrelevant or redundant with respect to the underlying cluster structure, thereby degrading clustering performance and interpretability. To address this, [38] proposed the 𝕊​ℝ​𝕌​𝕎 model, extending their earlier work [37], which assigns variables to one of four roles. Let 𝕊 denote the set of relevant clustering variables. Its complement, 𝕊c, comprises the irrelevant variables and is partitioned into subsets 𝕌 and 𝕎. Variables in 𝕌 are linearly explained by a subset ℝ⊆𝕊, while variables in 𝕎 are assumed independent of all relevant variables. This framework facilitates variable-specific interpretation and avoids over-penalization from complex covariance structures, as typically encountered in constrained GMMs [12]. The joint density under the 𝕊​ℝ​𝕌​𝕎 model, combining mixture components for clustering, regression, and independence, is defined as:

f𝕊​ℝ​𝕌​𝕎​(𝒚n∣K,m,r,l,𝕍,𝚯)=fclust​(𝒚n𝕊∣K,m,𝜶)​freg​(𝒚n𝕌∣r,𝒂+𝒚nℝ​𝜷,𝛀)​findep​(𝒚n𝕎∣l,𝜸,𝚪). (1)

Here, 𝕍=(𝕊,𝕌,ℝ,𝕎) denotes the variable partition, and 𝚯 is the full parameter set. The components are defined as: fclust:=fGMM, freg​(𝒚n𝕌∣r,𝒂+𝒚nℝ​𝜷,𝛀):=𝒩​(𝒚n𝕌∣𝒂+𝒚nℝ​𝜷,𝛀), with 𝒂∈ℝ1×|𝕌|, 𝜷∈ℝ|ℝ|×|𝕌|, and 𝛀∈ℝ|𝕌|×|𝕌| structured by r. The independent part is findep:=𝒩, with variance structure l and covariance 𝚪. Model selection proceeds by maximizing a BIC-type criterion:

critBIC​(K,m,r,l,𝕍)=BICclust​(𝒀𝕊∣K,m)+BICreg​(𝒀𝕌∣r,𝒀ℝ)+BICindep​(𝒀𝕎∣l), (2)

over (K,m,r,l,𝕍), where each term scores the corresponding model component. Although this approach offers fine-grained variable treatment, it can be computationally demanding in high dimensions due to stepwise selection for clustering and regression.

Penalized Log-Likelihood Methods for Simultaneous Clustering and Variable Selection. We follow the framework of [11], which extends [75, 13], introducing cluster-specific penalties via group-wise weighting matrices 𝑷k for adaptive regularization in Gaussian graphical mixture models. This method improves upon stepwise procedures [38, 37] by handling high-dimensional data more efficiently. The penalized log-likelihood is given by:

∑n=1Nln⁡[∑k=1Kπk​𝒩​(𝒚¯n∣𝝁k,𝚺k)]−λ​∑k=1K‖𝝁k‖1−ρ​∑k=1K∑d≠d′|(𝑷k⊙𝚺k−1)d​d′|. (3)

Here, 𝒚¯n denotes centered and scaled observations, and 𝑷k is typically chosen to be inversely proportional to initial partial correlations, enabling adaptive shrinkage akin to the adaptive lasso [77]. Given grids of regularization parameters 𝒢λ and 𝒢ρ, the EM algorithm from [75] is used to estimate the parameters 𝜶^​(λ,ρ)=(𝝅^​(λ,ρ),𝝁^1​(λ,ρ),…,𝝁^K​(λ,ρ),𝚺^1​(λ,ρ),…,𝚺^K​(λ,ρ)). For each variable d∈[D], clustering scores 𝒪K​(d) are computed, and the relevant set 𝕊 is formed from those that improve BIC, while 𝕌 captures non-informative variables. The final model is selected by maximizing critBIC​(K,m,r,l,𝕍(K,m,r)) in Equation 2, where 𝕍(K,m,r)=(𝕊(K,m),ℝ(K,m,r),𝕌(K,m),𝕎(K,m)).

3 Our Proposal: Variable Selection in Model-Based Clustering with MNAR

We propose a compact framework integrating adaptive variable selection with robust missing data handling in model-based clustering.

Global GMM Representation of 𝕊​ℝ​𝕌​𝕎 under MAR. A key aspect is expressing the observed-data likelihood in a tractable form. Under the MAR assumption, the 𝕊​ℝ​𝕌​𝕎 model admits a global GMM representation, extending the result of [39] for the 𝕊​ℝ model. Given the density in Equation 1 and Gaussian properties, the observed likelihood becomes

f​(𝒀o∣K,m,r,l,𝕍,𝚯)=∏n=1N(∑k=1Kπk​𝒩​(𝒚no∣𝝂~k,o,𝚫~k,o​o)), (4)

where 𝝂~k,o=(𝝂k,o𝜸o),𝚫~k,o​o=(𝚫k,o​o𝟎𝟎𝚪o​o). Here, 𝝂k,o,𝜸o and 𝚫k,o​o,𝚪o​o correspond to the block means and covariances from the 𝕊​ℝ model. The derivation extends the technical proof from [39] and is detailed in the supplement. This representation offers: (i) a unified GMM encoding the variable roles in 𝕊​ℝ​𝕌​𝕎 via observed-data parameters, (ii) compatibility with standard EM algorithms, enabling MAR-aware estimation and internal imputation that respects cluster structure, and (iii) theoretical guarantees of identifiability and consistency in selecting the true variable partition via Equation 2 (see Section 4).

Incorporating MNARz Mechanism into 𝕊​ℝ​𝕌​𝕎. To explicitly address MNAR data, we integrate the MNARz mechanism from [62], where missingness depends on the latent class membership. Estimation proceeds via EM, while identifiability and the transformation of MNARz into MAR on augmented data (i.e., original data with the missingness indicator matrix 𝑪) are retained (see Section 4). Let 𝒚n=(𝒚n𝕊,𝒚n𝕌,𝒚n𝕎) denote the full data for observation n, partitioned by role (Section 2), and let 𝒄n be its binary missingness pattern. With zn∈[K] denoting the latent class, the complete-data density under zn​k=1 is

f(𝒚n,𝒄n∣zn​k=1;𝜶k,𝝍k)=fclust(𝒚n𝕊∣𝜶k)freg(𝒚n𝕌∣𝒚nℝ;𝜽reg)×findep​(𝒚n𝕎∣𝜽indep)​fMNARz​(𝒄n∣zn​k=1;𝝍k),

where fMNARz​(𝒄n∣zn​k=1;𝝍k)=∏d=1Dρ​(𝝍k)cn​d​(1−ρ​(𝝍k))1−cn​d and ρ​(𝝍k) is the class-specific missingness probability. The model parameters are 𝚯=(π1,…,πK,{𝜶k}k=1K,𝝃={𝜽reg,𝜽indep},{𝝍k}k=1K). The observed-data log-likelihood is:

ℓ​(𝚯;𝒀,𝑪) =∑n=1Nlog​∑k=1Kπk​fko​(𝒚no;𝜶k,𝝃)​fc​(𝒄n;𝝍k), (5)

with the MNARz mechanism factored out due to independence. The complete-data log-likelihood is:

ℓcomp​(𝚯;𝒀,𝑪)=∑n=1N∑k=1Kzn​k​log⁡(πk​fk​(𝒚n;𝜶k,𝝃)​fc​(𝒄n;𝝍k)). (6)

To accommodate both MAR and MNAR patterns, we extend the model with a binary partition of variable indices: 𝒟MAR∪𝒟MNAR=[D],|𝒟MAR|=DM,|𝒟MNAR|=DM′. Variables in 𝒟MAR follow a MAR mechanism, and those in 𝒟MNAR follow MNARz:

fk​(𝒚n;𝜶k)=fclust​(𝒚n𝕊∣𝜶k)​freg​(𝒚n𝕌∣𝒚nℝ;𝜽reg)​findep​(𝒚n𝕎∣𝜽indep),𝝃=(𝜽reg,𝜽indep),
fcMAR​(𝒄n,M∣𝒚no;𝝍M)=∏d∈𝒟MARρd​(𝒚no;ψM​d)cn​d​(1−ρd​(𝒚no;ψM​d))1−cn​d,
fcMNARz​(𝒄n,N;𝝍k)=∏d∈𝒟MNARρ​(𝝍k)cn​d​(1−ρ​(𝝍k))1−cn​d.

The parameter blocks expand as: 𝚯=(π1,…,πK,{𝜶k}k=1K,𝝃,𝝍M,{𝝍k}k=1K). Let 𝒚n=(𝒚no,𝒚nm), and define 𝒄n,M=(cn​d)d∈𝒟MAR, 𝒄n,M′=(cn​d)d∈𝒟MNAR. Then, since the MNARz mask is independent of 𝒚n given k, the observed-data likelihood becomes:

ℓ​(𝚯;𝒀,𝑪)=∑n=1Nlog⁡[∑k=1Kπk​fcMNARz​(𝒄n,M′;𝝍k)​fk,Mo​(𝒚no;𝜶k,𝝍M)]. (7)

where fk,Mo​(𝒚no;𝜶k,𝝍M):=∫fk​(𝒚n;𝜶k)​fcMAR​(𝒄n,M∣𝒚no;𝝍M)​𝑑𝒚nm. The MNARz term factors out of the integral, preserving the ignorability property for that block [62]. Only the MAR component requires integration or imputation.

Adaptive Weighting for Penalty of Precision Matrix. To enhance performance, we replace the inverse partial correlation scheme with a spectral-based computation of π𝐤 in Equation 3. Starting from the initial precision matrix 𝚿^k(0), we construct an unweighted, undirected graph 𝒢(k)=(𝒱,ℰ(k)), where 𝒱=[D], and include edge (i,j) in ℰ(k) if ∥𝝍^k,i​j(0)∥>𝚪a​d​j. Let 𝐀(k) and 𝐃(k) denote the adjacency and degree matrices, respectively. The symmetrically normalized Laplacian is then defined as 𝐋s​y​m(k)=𝐈−(𝐃(k))−1/2​𝐀(k)​(𝐃(k))−1/2, which offers a scale-invariant measure of connectivity. We aim to shrink 𝚿k toward a diagonal target (i.e., an empty graph), with 𝐋target=𝟎. The spectral distance between 𝚿^k(0) and this target is DLS​(𝚿k(0))=∥spec​(𝐋symk)∥2, and the corresponding adaptive weights are given by 𝐏k,i​j=(DLS​(𝚿k(0))+ϵ)−1.

Parameter Estimation. Following [75], we adopt similar parameter estimation for variable ranking and focus here on estimating 𝕊​ℝ​𝕌​𝕎 under the MNARz mechanism using the EM algorithm. Other estimation details are provided in the Supplement. Let each observation be partitioned as 𝒚n=(𝒚no,𝒚nm), and define 𝝃=(𝜽reg,𝜽indep). Since fc is independent of 𝒚n, the observed-data log-likelihood becomes:

ℓ​(𝚯;𝒀o,𝑪)=∑n=1Nlog​∑k=1Kπk​fc​(𝒄n;𝝍k)​fko​(𝒚no;𝜶k,𝝃),

where fko​(𝒚no;𝜶k,𝝃):=∫fclust​(𝒚nS∣𝜶k)​freg​(𝒚nU∣𝒚nR;𝜽reg)​findep​(𝒚nW∣𝜽indep)​𝑑𝒚nm.

The E-step computes responsibilities:

tn​k(t)=πk(t−1)​fko​(𝒚no;𝜶k(t−1),𝝃(t−1))​fc​(𝒄n;𝝍k(t−1))∑l=1Kπl(t−1)​flo​(𝒚no;𝜶l(t−1),𝝃(t−1))​fc​(𝒄n;𝝍l(t−1)).

The Q-function to be maximized in the M-step is:

Q​(𝚯;𝚯(t−1)) =∑n=1N∑k=1Ktn​k(t)​{log⁡πk+gy​(𝜶k)+gc​(𝝍k)}.
gy​(𝜶k) =𝔼𝚯(t−1)​[log⁡fk​(𝒚n;𝜶k,𝝃k)∣𝒚no,𝒄n,zn​k=1],
gc​(𝝍k) =𝔼𝚯(t−1)​[log⁡fc​(𝒄n∣𝒚n;𝝍k)∣𝒚no,𝒄n,zn​k=1].

In the M-step, we update:

πk(t) =1N​∑ntn​k(t),𝜶k(t)=arg⁡max𝜶k​∑ntn​k(t)​gy​(𝜶k),
𝝍k(t) =arg⁡max𝝍k​∑ntn​k(t)​gc​(𝝍k),𝝃(t)=arg⁡max𝝃​∑n,ktn​k(t)​gy​(𝝃).

We now summarize the workflow; the penalty acts only in Stage A (ranking), while Stage B performs unpenalized SRUW model selection by a single pass over the ranked list.

Algorithm 1 High level two-stage SRUW-MNARz procedure

Input: incomplete data 𝒀, mask 𝑪, model grid (K,m); regularization grids 𝒢λ (and 𝒢ρ).
Stage A: Ranking (penalized GMM on a fast-imputed data)

  1. 1.

    Produce a fast single imputation 𝒀~ only to enable ranking.

  2. 2.

    For each (λ,ρ)∈𝒢λ×𝒢ρ, fit the penalized GMM to 𝒀¯ where 𝒀¯ denotes the centered/scaled 𝒀~; record μ^k​d​(λ,ρ).

  3. 3.

    Compute ranking score 𝒪K​(d) for each variable d by counting along the path how often {μ^k​d​(λ,ρ)}k remain nonzero; sort variables by 𝒪K​(d) (see [13] for details).

Stage B: Role assignment (SRUW on original incomplete data)

  1. 1.

    Traverse ranked variables once (similar to [13]); at each step, fit unpenalized SRUW-MNARz on incomplete 𝒀 and decide the role (S/U/W) of the new variable using BIC criterion using Equation 2 as in [38]

  2. 2.

    Return (𝕊^,ℝ^,𝕌^,𝕎^), K^,m^, and parameter estimates.

Output: final role sets and estimates.

Practical Initialization. Like all mixture EM algorithms, ours is sensitive to initialization. We use a robust warm start: (i) fast single imputation 𝒀~ (only for ranking); (ii) K-means++/hierarchical clustering/random start on 𝒀~ for (πk,𝝁k,𝚺k); (iii) given initial partitions from step (ii) 𝚿k(0) set to diagonal estimates using samples assigned to each cluster (cluster-aware estimates); (iv) 𝝍k initialized from class-wise missing rates after a soft E-step. We then run a λ-path (as illustrated below) for Stage A, followed by unpenalized SRUW estimation for Stage B. Per-K penalty grids are constructed by data-dependent upper bounds and a geometric path: For λ, with the hard partition 𝒁(0) from initialization,

λmax=maxk∈[K],d∈[D]⁡|(𝒁(0)⊤​𝑿)k​d|,λℓ=λmax​(λminλmax)ℓ−1L−1,λmin=ξ​λmax,ℓ=[L].

This matches the KKT threshold at which the ℓ1 penalty zeros all 𝝁k entries. For ρ, let 𝑺k be the (hard-label) empirical covariance and note our Stage-A glasso objective uses the coefficient (2​ρ/nk) in front of the weighted ℓ1 term. A diagonal-forcing threshold is

ρmax=maxk⁡maxi≠j⁡nk​|(𝑺k)i​j|(𝑷k)i​j,

followed by the same L-point geometric path ρℓ=ρmax​(ρmin/ρmax)(ℓ−1)/(L−1) with ρmin=ξ​ρmax. We set diag⁡(𝑷k)=𝟎 when we do not impose a penalty on diagonal entries.

4 Theoretical Properties

In this section, we present theoretical guarantees for our extended framework, including identifiability and selection consistency under the MNARz mechanism, and extend the results to the general 𝕊​ℝ​𝕌​𝕎 model. We also show that these guarantees apply to the regularized approach of [13]. Given the known issues in finite mixture models, such as degeneracy and non-identifiability, we adopt standard assumptions, largely consistent with those in [38]. For any configuration (𝕊,ℝ,𝕌,𝕎), let 𝜽(𝕊,ℝ,𝕌,𝕎)∗ denote the true parameter and 𝜽^(𝕊,ℝ,𝕌,𝕎) its maximum likelihood estimator.

Theorem 1 (Informal: Identifiability of the 𝕊​ℝ​𝕌​𝕎 Model).

Let (K,m,r,l,𝕍) and (K∗,m∗,r∗,l∗,𝕍∗) denote two models under the MNARz mechanism. Let 𝚯(K,m,r,l,𝕍)⊆𝚼(K,m,r,l,𝕍) denote the parameter space such that each element 𝛉=(𝛂,𝐚,𝛃,𝛀,𝛄,𝚪) satisfies the following: (i) the component-specific parameters (𝛍k,𝚺k) are distinct and satisfy, for all s⊆𝕊, there exist 1≤k<k′≤K such that 𝛍k,s¯|s≠𝛍k′,s¯|sor𝚺k,s¯|s≠𝚺k′,s¯|sor𝚺k,s¯​s¯|s≠𝚺k′,s¯​s¯|s, where s¯ denotes the complement of s in 𝕊; (ii) if 𝕌≠∅, then for every j∈ℝ, there exists u∈𝕌 such that 𝛃u​j≠0, and each u∈𝕌 affects at least one j∈ℝ with 𝛃u​j≠0; and (iii) the parameters 𝚪 and 𝛄 exactly respect the structural forms r and l, respectively. If there exist 𝛉∈𝚯(K,m,r,l,𝕍) and 𝛉∗∈𝚯(K∗,m∗,r∗,l∗,𝕍∗) such that f(⋅∣K,m,r,l,𝕍,𝛉)=f(⋅∣K∗,m∗,r∗,l∗,𝕍∗,𝛉∗), then (K,m,r,l,𝕍)=(K∗,m∗,r∗,l∗,𝕍∗) and 𝛉=𝛉∗.

Theorem 2 (Informal: BIC consistency for the 𝕊​ℝ​𝕌​𝕎 model).

Assume the MNARz mechanism, the regularity conditions specified in the Supplementary Material, and the existence of a unique tuple (K0,m0,r0,l0,𝕍0) such that the data-generating distribution satisfies h=f(⋅∣𝛉(K0,m0,r0,l0,𝕍0)∗) for some true parameter 𝛉∗ where h is the density function of the sample 𝐘 and the model specification (K0,m0,r0,l0,𝕍0) is assumed to be known. Let (K0,m0,r0,l0) be fixed. Then the variable selection procedure that selects the subsets (𝕊^,ℝ^,𝕌^,𝕎^) by maximizing the BIC criterion is consistent, in the sense that ℙ​((𝕊^,ℝ^,𝕌^,𝕎^)=(𝕊0,ℝ0,𝕌0,𝕎0))→N→∞1.

Theorem 3 (Informal: Selection consistency of the two-stage procedure).

Under certain regularity conditions in Supplementary Material, the backward stepwise selection procedure in the 𝕊​ℝ​𝕌​𝕎 model, guided by the variable ranking stage, consistently recovers the true relevant variable set 𝕊0∗ with high probability. That is, ℙ​(𝕊^∗=𝕊0∗)→N→∞1.

Proof sketch. The proof of Theorem 3 proceeds in three main steps: 1) We first establish the consistency of the penalized M-estimator used in the variable ranking stage, specifically the penalized GMMs. 2) This consistency implies that the variables are ranked in a manner that, with high probability, separates relevant variables from noise variables. 3) We then show that the subsequent BIC-based role determination step, when applied to this consistent ranking, correctly recovers the true variable roles. A key component of the proof involves establishing consistency in parameter estimation. In particular, we derive a Restricted Strong Convexity (RSC)-type condition for the negative GMM log-likelihood. More precisely, let ℒN​(𝜶) denote the empirical objective function based on sample size N, and let 𝜶⋆ be the true parameter. There exist constants γ1,γ2>0 such that for every perturbation 𝚫 in a restricted neighborhood of 𝜶⋆, ℒN​(𝜶⋆+𝚫)−ℒN​(𝜶⋆)−⟨∇ℒN​(𝜶⋆),𝚫⟩≥γ12​∥𝚫𝝁∥22+γ22​∥𝚫𝚺−1∥F2−τN​(𝚫), where τN​(𝚫) is a tolerance term that decays at rate N−1/2. A complete and rigorous proof is provided in the Supplementary Material.

Implications of the theorems. In the framework of [13], the parameter c in the BIC-based selection step (representing the number of consecutive non-positive BIC differences allowed before terminating inclusion into 𝕊 or 𝕎) governs the balance between false negatives (omitting true variables) and false positives (including irrelevant variables). Theorem 3 implies that for finite samples, choosing a moderately large c, such as c≈log⁡(D) (e.g., 3 to 5, as commonly practiced), provides robustness against spurious BIC fluctuations caused by noise variables. In particular, if the gap in the 𝒪K​(d) scores between the last true signal and the first noise variable is sufficiently pronounced, the exact value of c is not overly critical for ensuring asymptotic consistency, as long as it adequately accounts for statistical noise in the BIC differences sequence. However, an excessively large c could lead to the inclusion of noise variables, especially if the variable ranking is imperfect.

5 Experiments

Our primary objective in this section is to evaluate the clustering quality, imputation accuracy, as well as variable-role recovery under incomplete data scenarios. For the simulated dataset, we compare our method with its direct predecessors in [37, 13], denoting them as Clustvarsel and Selvar respectively, and design each of them to be a pre-imputation+clustering pipeline with the imputation method to be Random Forest [10] using missRanger package [40, 63]. For [37], we use the algorithm in [59] with forward direction and headlong process, with rationale being explained in [13].

We also benchmark with VarSelLCM illustrated in [35] since the model itself also recasts the variable selection problem into a model-selection one. For this case study, we adopt the same real high-dimensional dataset on transcriptome as in [37, 39, 13] and interpret the result in comparison with previous findings in [37, 39]. To measure imputation accuracy, we use WNRMSE - a weighted version of NRMSE where the weights are deduced from the proportion of missingness in each group, while clustering performance is measured by ARI. The variable selection will be assessed by the frequency the model chooses the correct set of relevant variables, which are controlled in the simulation.

Simulated Dataset. We generate two synthetic data sets with sample size n=2000 for variable/model selection with missing values treatment benchmarking. The first data set mimics [39] where each observation 𝒚n∈ℝ7 is drawn from a four-component Gaussian mixture on the clustering block 𝕊={𝒚1,𝒚2,𝒚3}. The components are equiprobable and have means (0,0,0)′,(−6,6,0)′,(0,0,6)′,(−6,6,6)′ with a diagonal covariance matrix 𝚺=diag​(6​2,1,2)​σscale2. The redundant block 𝕌={𝒚4,𝒚5} is generated by a linear regression 𝕌=(−1,2)′+𝐗S​[(0.5,2)′,(1,0)′]+ϵ, where ϵ∼𝒩​(𝟎,𝐑​(π/𝟔)​diag​(𝟏,𝟑)​𝐑​(π/𝟔)′) and 𝐑​(𝜽) denotes the usual planar rotation. The remaining W variables are iid 𝒩​(0,1) noise. In the second design, an observation 𝒚n∈ℝ14. The clustering block now consists of (𝒚1,𝒚2), drawn from an equiprobable four-component Gaussian mixture with means (0,0),(4,0),(0,2),(4,2) and common covariance 0.5​𝐈𝟐. Conditional on (𝒚1,𝒚2), the vector 𝒚3:14 follows a linear model 𝒚3:14=𝜶+(𝒚1,𝒚2)​𝜷+𝜺, where 𝜶, the 2×12 coefficient matrix 𝜷 and the diagonal covariance 𝛀 are varied across eight sub-scenarios: scenario 1 contains only pure noise, scenarios 2​-​7 introduce more regressors, and scenario 8 adds intercept shifts together with three additional pure-noise variables 𝒚12:14. In our setting, we consider K∈{2,3,4} and every method except VarSelLCM is allowed to learn its own covariance structure. We vary the missing ratio {0.05,0.1,0.2,0.3,0.5} under both MAR and MNAR patterns. For the second simulated dataset, we temporarily consider only scenario 8. Moreover, for Selvar and SelvarMNARz, we fix c=2 and compute 𝑷k with the spectral distance. Further details of the synthetic-data generation scheme are provided in the supplement. From Figure 1, we observe that SelvarMNARz delivers the highest ARI and the lowest WNRMSE across every missing-data level in both scenarios. Its performance declines only modestly as missingness increases, while the impute-then-cluster baselines deteriorate sharply, especially at 30% and 50% missingness under MNAR, pointing to the benefit of modeling the missing-data mechanism within the clustering algorithm. To statistically validate that our model exhibits a significantly slower performance decline under high missingness, we conducted one-sided Welch’s t-tests (Bonferroni-corrected) at the 50% MNAR level. As shown in Table 1, the ARI of SelvarMNARz is significantly greater than all baseline models (p<0.001) across 20 replications.

Table 1: Welch’s t-test for ARI (MNAR, 50% missingness, α=0.05).
Model Mean ARI Std Comparison p-value Signif.
SelvarMNARz 0.511 0.052 – – NA
Clustvarsel 0.363 0.088 vs. SelvarMNARz <0.001 ∗⁣∗⁣∗⁣∗
Selvar 0.348 0.108 vs. SelvarMNARz <0.001 ∗⁣∗⁣∗⁣∗
VarSelLCM 0.344 0.101 vs. SelvarMNARz <0.001 ∗⁣∗⁣∗⁣∗

While all methods show a monotonic performance drop as missingness increases, these results confirm that the superior performance of SelvarMNARz is not a random artifact but a statistically significant improvement, demonstrating its robustness in the challenging missing data scenarios. Figure 2 shows that SelvarMNARz consistently recovers the true cluster component and the correct set of clustering variables, regardless of the missing-data mechanism or rate. On Dataset 2, every method succeeds under both MAR and MNAR, whereas on Dataset 1, the performance of Clustvarsel and VarSelLCM drops as missingness increases, leading to the omission of important variables. The contrast reflects the data-generating designs: Dataset 2 includes pure-noise variables 𝕎 within the 𝕊​ℝ​𝕌​𝕎 framework, whereas Dataset 1 aligns more closely with the assumptions behind Clustvarsel and VarSelLCM. Finally, workflows that rely on external imputation misidentify relevant variables once missingness exceeds 20% under either missing-data scenarios.

Refer to caption
Refer to caption
Figure 1: Comparison of four models under MAR and MNAR mechanisms over 20 replications; for ARI/WNRMSE, higher/lower boxplots indicate better performance.
Refer to caption
Figure 2: Proportions choosing correct relevant variables and cluster components over 20 replications.
Refer to caption
Figure 3: Mean expression profiles across 18 clusters. Light region indicates irrelevant P.
Table 2: Size of each mixture component and R2 obtained by regressing the 𝕌-block on the 𝕊-block within that cluster. Values extremely close to 1 (marked by ∗) occur because the 𝕌-block is nearly constant in those groups.
Cluster # Genes 𝑹𝟐
1 715 0.02
2 11 1.00∗
3 6 1.00∗
4 93 0.45
5 124 0.12
6 9 1.00∗
Cluster # Genes 𝑹𝟐
7 36 0.68
8 13 1.00∗
9 15 1.00∗
10 29 0.59
11 71 0.43
12 24 0.91
Cluster # Genes 𝑹𝟐
13 46 0.57
14 19 1.00∗
15 6 1.00∗
16 16 1.00∗
17 24 0.90
18 10 1.00∗

Transcriptome Data. We re-examine the 1,267-gene Arabidopsis thaliana transcriptome dataset previously analyzed using SelvarClust and SelvarClustMV in [37, 39] to illustrate clustering with variable selection. Full experimental details are provided in the Supplementary Material; here, we summarize the settings and key findings. We fitted SelvarMNARz for K∈{2,…,20} clusters with c=5, spectral distance weights 𝑷k, and a πk​L​C mixture structure as in [39]. For conciseness, we refer to Project 1 as P1, Project 2 as P2, and so on. As shown in Figure 3, the proposed framework partitions the transcriptome into 18 clusters and reaffirms P1-P4 as core drivers, while globally assigning P5-P7 as being explained by these drivers. This finding is consistent with [39], which also identifies P2 as relevant for clustering, supporting the hypothesis that iron signaling responses are critical for defining co-expression groups when a broader gene set (including entries with missing data) is considered. Additionally, all three methods agree that P5 is not a primary clustering factor and is more likely a correlated effect explained by broader cellular processes. However, a key departure lies in the reclassification of P6 and P7, previously identified as informative clustering drivers [37, 39], into a redundant 𝕌 block, explained by P1-P4. This suggests that, after accounting for MNAR effects, variation in P6-P7 is largely predictable from the core axes P1-P4, and P5-P7 all measure late-stage stress outputs. Further analysis of cluster sizes and the per-cluster R2 values, shown in Table 2, validates the selected variables when regressing P5-P7 on P1-P4 within each cluster. Six biologically coherent clusters (6, 7, 8, 10, 12, 18) exhibit R2>0.60, indicating that late-stress variation is well encoded by P1-P4. Cluster 1, although large, is transcriptionally flat across both 𝕊 and 𝕌; its low R2=0.02 thus reflects a vanishing denominator rather than hidden structure. True residual structure (unexplained 𝕌 signal) is concentrated in two groups: clusters 5 and 11, which are exactly where the 18th mixture component forms.

6 Conclusion, Limitations and Perspectives

We proposed a unified model-based clustering framework that jointly performs variable selection and handles MNAR data by combining adaptive penalization with explicit modeling of missingness-latent class dependencies. This enhances clustering flexibility and robustness in high-dimensional settings like transcriptomics. The method enjoys theoretical guarantees, including selection and asymptotic consistency, and shows strong empirical performance. However, it is currently limited to continuous data via Gaussian mixture models, making it less suitable for categorical or mixed-type variables. Extending the framework to categorical settings is a natural next step. For instance, Dean and Raftery [15] introduced variable selection for latent class models, later refined by Fop et al. [20] using a statistically grounded stepwise approach. Bontemps and Toussile [8] proposed a mixture of multivariate multinomial distributions with combinatorial variable selection and slope heuristics for penalty calibration. Building upon these, a future extension of our framework could incorporate discrete data distributions, adapt penalization accordingly, and model missingness mechanisms using latent class structures. Further directions include mixed-type data extensions and scalability improvements to accommodate larger, more heterogeneous datasets.

Acknowledgments

TrungTin Nguyen and Christopher Drovandi acknowledge funding from the Australian Research Council Centre of Excellence for the Mathematical Analysis of Cellular Systems (CE230100001). Christopher Drovandi was also supported by an Australian Research Council Future Fellowship (FT210100260). This paper was also supported by AISIA Research Lab, Vietnam. The authors would like to express their special thanks to the authors of [37, 39], in particular Professor Cathy Maugis-Rabusseau and Dr. Marie-Laure Martin, for providing the Arabidopsis thaliana transcriptome dataset.

References

  • [1] H. Akaike. A new look at the statistical model identification. IEEE transactions on automatic control, 19(6):716–723, 2003. ↩ 1 G.2
  • [2] E. S. Allman, C. Matias, and J. A. Rhodes. Identifiability of parameters in latent structure models with many observed variables. Annals of Statistics, 37(6A):3099–3132, 2009. ↩ D.1
  • [3] J. L. Andrews and P. D. McNicholas. Variable selection for clustering and classification. Journal of Classification, 31(2):136–153, 2014. ↩ 1 G.2
  • [4] I. Balelli, A. Sportisse, F. Cremonesi, P.-A. Mattei, and M. Lorenzi. Fed-miwae: Federated imputation of incomplete data via deep generative models. arXiv preprint arXiv:2304.08054, 2023. ↩ 1 G.3
  • [5] E. Barut, J. Fan, and A. Verhasselt. Conditional sure independence screening. Journal of the American Statistical Association, 111(515):1266–1277, 2016. ↩ 1
  • [6] J.-P. Baudry, C. Maugis, and B. Michel. Slope heuristics: overview and implementation. Statistics and Computing, 22(2):455–470, 2012. ↩ 2
  • [7] C. Biernacki, G. Celeux, and G. Govaert. Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(7):719–725, 2000. ↩ 2
  • [8] D. Bontemps and W. Toussile. Clustering and variable selection for categorical multivariate data. Electronic Journal of Statistics, 7:2344–2371, Jan. 2013. ↩ 6 G.2
  • [9] L. Breiman. Better subset regression using the nonnegative garrote. Technometrics, 37(4):373–384, 1995. ↩ 1 G.2
  • [10] L. Breiman. Random forests. Machine Learning, 45:5–32, 2001. ↩ 5
  • [11] A. Casa, A. Cappozzo, and M. Fop. Group-wise shrinkage estimation in penalized model-based clustering. Journal of Classification, 39(3):648–674, 2022. ↩ 2 A
  • [12] G. Celeux and G. Govaert. Gaussian parsimonious clustering models. Pattern Recognition, 28(5):781–793, 1995. Publisher: Elsevier. ↩ 2
  • [13] G. Celeux, C. Maugis-Rabusseau, and M. Sedki. Variable selection in model-based clustering and discriminant analysis with a regularization approach. Advances in Data Analysis and Classification, 13:259–278, 2019. ↩ 1 2 3 4 5 A F.3 F.5 G.2
  • [14] M. C. Chong, H. D. Nguyen, and TrungTin Nguyen. Risk Bounds for Mixture Density Estimation on Compact Domains via the h-Lifted Kullback–Leibler Divergence. Transactions on Machine Learning Research, 2024. ↩ 2
  • [15] N. Dean and A. E. Raftery. Latent class analysis variable selection. Annals of the Institute of Statistical Mathematics, 62:11–35, 2010. ↩ 6 G.2
  • [16] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1–38, 1977. ↩ 1 G.1
  • [17] D. Do, L. Do, S. A. McKinley, J. Terhorst, and X. Nguyen. Dendrogram of mixing measures: Learning latent hierarchy and model selection for finite mixture models. arXiv preprint arXiv:2403.01684, 2024. ↩ 2
  • [18] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression. Annals of Statistics, pages 407–451, 2004. ↩ 1 G.2
  • [19] J. Fan and J. Lv. Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society Series B: Statistical Methodology, 70(5):849–911, 2008. ↩ 1
  • [20] M. Fop, K. M. Smart, and T. B. Murphy. Variable selection for latent class analysis with application to low back pain diagnosis. The Annals of Applied Statistics, pages 2080–2110, 2017. ↩ 6 G.2
  • [21] F. Forbes, H. D. Nguyen, T. Nguyen, and J. Arbel. Mixture of expert posterior surrogates for approximate Bayesian computation. In JDS 2022 - 53èmes Journées de Statistique de la Société Française de Statistique (SFdS), Lyon, France, June 2022. ↩ 2
  • [22] F. Forbes, H. D. Nguyen, T. Nguyen, and J. Arbel. Summary statistics and discrepancy measures for approximate Bayesian computation via surrogate posteriors. Statistics and Computing, 32(5):85, Oct. 2022. ↩ 2
  • [23] R. Foygel and M. Drton. Extended Bayesian Information Criteria for Gaussian Graphical Models. In J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and A. Culotta, editors, Advances in Neural Information Processing Systems, volume 23. Curran Associates, Inc., 2010. ↩ 2
  • [24] S. Frühwirth-Schnatter, C. Pamminger, A. Weber, and R. Winter-Ebmer. Labor market entry and earnings dynamics: Bayesian inference using mixtures-of-experts Markov chain clustering. Journal of Applied Econometrics, 27(7):1116–1137, 2012. ↩ 2
  • [25] C. R. Genovese and L. Wasserman. Rates of convergence for the Gaussian mixture sieve. The Annals of Statistics, 28(4):1105 – 1127, 2000. ↩ 2
  • [26] N. Ho and X. Nguyen. Convergence rates of parameter estimation for some weakly identifiable finite mixtures. The Annals of Statistics, 44(6):2726 – 2755, 2016. ↩ 2
  • [27] N. Ho and X. Nguyen. On strong identifiability and convergence rates of parameter estimation in finite mixtures. Electronic Journal of Statistics, 10(1):271–307, 2016. ↩ 2
  • [28] A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1):55–67, 1970. ↩ 1 G.2
  • [29] A. Jasra, C. C. Holmes, and D. A. Stephens. Markov Chain Monte Carlo Methods and the Label Switching Problem in Bayesian Mixture Modeling. Statistical Science, 20(1):50 – 67, 2005. ↩ 1 G.1
  • [30] S. Kim, M. G. Tadesse, and M. Vannucci. Variable selection in clustering via dirichlet process mixture models. Biometrika, 93(4):877–893, 2006. ↩ 1
  • [31] J. Kwon and C. Caramanis. The em algorithm gives sample-optimality for learning mixtures of well-separated gaussians. In Conference on Learning Theory, pages 2425–2487. PMLR, 2020. ↩ F.4
  • [32] M. H. Law, M. A. Figueiredo, and A. K. Jain. Simultaneous feature selection and clustering using mixture models. IEEE transactions on pattern analysis and machine intelligence, 26(9):1154–1166, 2004. ↩ 1 G.2
  • [33] R. J. Little and D. B. Rubin. Statistical analysis with missing data. John Wiley & Sons, 2019. ↩ A
  • [34] P.-L. Loh and M. J. Wainwright. Regularized M-estimators with Nonconvexity: Statistical and Algorithmic Theory for Local Optima. Journal of Machine Learning Research, 16(19):559–616, 2015. ↩ A
  • [35] M. Marbac and M. Sedki. Variable selection for mixed data clustering: A model-based approach. arXiv preprint arXiv:1703.02293, 2017. ↩ 5 A
  • [36] P.-A. Mattei and J. Frellsen. MIWAE: Deep Generative Modelling and Imputation of Incomplete Data Sets. In K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 4413–4423. PMLR, June 2019. ↩ 1 G.3
  • [37] C. Maugis, G. Celeux, and M.-L. Martin-Magniette. Variable selection for clustering with gaussian mixture models. Biometrics, 65(3):701–709, 2009. ↩ 2 5 ↑ B D.1 F.6
  • [38] C. Maugis, G. Celeux, and M.-L. Martin-Magniette. Variable selection in model-based clustering: A general variable role modeling. Computational Statistics & Data Analysis, 53(11):3872–3882, 2009. ↩ 1 2 3 4 A B D.1 D.2 G.2
  • [39] C. Maugis-Rabusseau, M.-L. Martin-Magniette, and S. Pelletier. Selvarclustmv: Variable selection approach in model-based clustering allowing for missing values. Journal de la Société Française de Statistique, 153(2):21–36, 2012. ↩ 1 2 3 5 ↑ F.6
  • [40] M. Mayer. missRanger: Fast Imputation of Missing Values, 2025. R package version 2.6.2, https://mayer79.github.io/missRanger/. ↩ 5
  • [41] R. Mazumder and T. Hastie. Exact covariance thresholding into connected components for large-scale graphical lasso. Journal of Machine Learning Research, 13(27):781–794, 2012. ↩ F.4
  • [42] G. J. McLachlan and K. E. Basford. Mixture Models: Inference and Applications to Clustering, volume 38. Marcel Dekker, 1988. ↩ 1
  • [43] G. J. McLachlan and D. Peel. Finite mixture models. John Wiley & Sons, 2004. ↩ 1
  • [44] S. Negahban and M. J. Wainwright. Restricted strong convexity and weighted matrix completion: Optimal bounds with noise. Journal of Machine Learning Research, 13(53):1665–1697, 2012. ↩ A D.3
  • [45] D. N. Nguyen and Z. Li. Joint learning of Gaussian graphical models in heterogeneous dependencies of high-dimensional transcriptomic data. In The 16th Asian Conference on Machine Learning (Conference Track), 2024. ↩ 2
  • [46] H. Nguyen, T. Nguyen, and N. Ho. Demystifying Softmax Gating in Gaussian Mixture of Experts. In Advances in Neural Information Processing Systems, Dec. 2023. ↩ 2
  • [47] H. Nguyen, T. Nguyen, K. Nguyen, and N. Ho. Towards Convergence Rates for Parameter Estimation in Gaussian-gated Mixture of Experts. In S. Dasgupta, S. Mandt, and Y. Li, editors, Proceedings of The 27th International Conference on Artificial Intelligence and Statistics, volume 238 of Proceedings of Machine Learning Research, pages 2683–2691. PMLR, May 2024. ↩ 2
  • [48] H. D. Nguyen, T. Nguyen, and F. Forbes. Bayesian Likelihood Free Inference using Mixtures of Experts. In 2024 International Joint Conference on Neural Networks (IJCNN), pages 1–8, 2024. ↩ 2
  • [49] T. Nguyen, F. Chamroukhi, H. D. Nguyen, and G. J. McLachlan. Approximation of probability density functions via location-scale finite mixtures in Lebesgue spaces. Communications in Statistics - Theory and Methods, 52(14):5048–5059, 2023. ↩ 2
  • [50] T. Nguyen, D. N. Nguyen, H. D. Nguyen, and F. Chamroukhi. A non-asymptotic theory for model selection in high-dimensional mixture of experts via joint rank and variable selection. In AJCAI Australasian Joint Conference on Artificial Intelligence 2023, Dec. 2023. ↩ 2
  • [51] T. Nguyen, H. D. Nguyen, F. Chamroukhi, and F. Forbes. A non-asymptotic approach for model selection via penalization in high-dimensional mixture of experts models. Electronic Journal of Statistics, 16(2):4742 – 4822, 2022. ↩ 2
  • [52] T. Nguyen, H. D. Nguyen, F. Chamroukhi, and G. J. McLachlan. Approximation by finite mixtures of continuous density functions that vanish at infinity. Cogent Mathematics & Statistics, 7(1):1750861, Jan. 2020. ↩ 2
  • [53] X. Nguyen. Convergence of latent mixing measures in finite and infinite mixture models. The Annals of Statistics, 41(1):370–400, 2013. ↩ 2
  • [54] A. E. Raftery and N. Dean. Variable selection for model-based clustering. Journal of the American Statistical Association, 101(473):168–178, 2006. ↩ 1 G.2
  • [55] A. Rakhlin, D. Panchenko, and S. Mukherjee. Risk bounds for mixture density estimation. ESAIM: PS, 9:220–229, 2005. ↩ 2
  • [56] D. B. Rubin. Inference and missing data. Biometrika, 63(3):581–592, 1976. ↩ A
  • [57] D. B. Rubin. An overview of multiple imputation. In Proceedings of the survey research methods section of the American statistical association, volume 79, page 84. Citeseer, 1988. ↩ 1 G.3
  • [58] G. Schwarz. Estimating the dimension of a model. The annals of statistics, pages 461–464, 1978. ↩ 1 2 G.2
  • [59] L. Scrucca and A. E. Raftery. clustvarsel: A package implementing variable selection for gaussian model-based clustering in r. Journal of Statistical Software, 84:1–28, 2018. ↩ 1 5 G.2
  • [60] W. Shen, S. T. Tokdar, and S. Ghosal. Adaptive Bayesian multivariate density estimation with Dirichlet mixtures. Biometrika, 100(3):623–640, 2013. ↩ 2
  • [61] C.-Y. Sin and H. White. Information criteria for selecting possibly misspecified parametric models. Journal of Econometrics, 71(1):207–225, 1996. ↩ 2
  • [62] A. Sportisse, M. Marbac, F. Laporte, G. Celeux, C. Boyer, J. Josse, and C. Biernacki. Model-based clustering with missing not at random data. Statistics and Computing, 34(4):135, 2024. ↩ 1 3 A F.2 F.5 G.3
  • [63] D. J. Stekhoven and P. Bühlmann. Missforest—non-parametric missing value imputation for mixed-type data. Bioinformatics, 28(1):112–118, 2012. ↩ 1 5 G.3
  • [64] M. Stephens. Bayesian analysis of mixture models with an unknown number of components—an alternative to reversible jump methods. The Annals of Statistics, 28(1):40 – 74, 2000. ↩ 1
  • [65] M. Stephens. Dealing with label switching in mixture models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(4):795–809, 2000. ↩ G.1
  • [66] N. Städler, P. Bühlmann, and S. van de Geer. l1-penalization for mixture regression models. TEST, 19(2):209–256, Aug. 2010. ↩ A
  • [67] M. G. Tadesse, N. Sha, and M. Vannucci. Bayesian variable selection in clustering high-dimensional data. Journal of the American Statistical Association, 100(470):602–617, 2005. ↩ 1
  • [68] H. Teicher. Identifiability of finite mixtures. The Annals of Mathematical Statistics, 34(4):1265–1269, December 1963. ↩ D.1
  • [69] T. Thai, T. Nguyen, D. Do, N. Ho, and C. Drovandi. Model Selection for Gaussian-gated Gaussian Mixture of Experts Using Dendrograms of Mixing Measures. arXiv preprint arXiv:2505.13052, 2025. ↩ 2
  • [70] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267–288, 1996. ↩ 1 A G.2
  • [71] S. Van Buuren and K. Groothuis-Oudshoorn. mice: Multivariate imputation by chained equations in r. Journal of statistical software, 45:1–67, 2011. ↩ 1 G.3
  • [72] J. Westerhout, T. Nguyen, X. Guo, and H. D. Nguyen. On the Asymptotic Distribution of the Minimum Empirical Risk. In Forty-first International Conference on Machine Learning, 2024. ↩ 2
  • [73] D. M. Witten, J. H. Friedman, and N. Simon. New insights and faster computations for the graphical lasso. Journal of Computational and Graphical Statistics, 20(4):892–900, 2011. ↩ F.4
  • [74] Y. Zhao and Q. Long. Multiple imputation in the presence of high-dimensional data. Statistical Methods in Medical Research, 25(5):2021–2035, 2016. ↩ 1 G.3
  • [75] H. Zhou, W. Pan, and X. Shen. Penalized model-based clustering with unconstrained covariance matrices. Electronic Journal of Statistics, 3(none):1473 – 1496, 2009. Publisher: Institute of Mathematical Statistics and Bernoulli Society. ↩ 1 2 3 A G.2
  • [76] M. Zhou et al. Global convergence of gradient em for over-parameterized gaussian mixtures. arXiv preprint arXiv:2506.06584, 2025. ↩ F.4
  • [77] H. Zou. The adaptive lasso and its oracle properties. Journal of the American Statistical Association, 101(476):1418–1429, 2006. ↩ 2

Supplementary Materials for
“A Unified Framework for Variable Selection in Model-Based Clustering with Missing Not at Random”

In this supplementary material, we first provide additional details on the main challenges as well as the technical and computational contributions in Appendix A, with the aim of enhancing the reader’s understanding of our key theoretical and methodological developments. Next, we summarize the definitions and theoretical results for the 𝕊​ℝ and 𝕊​ℝ​𝕌​𝕎 models in Appendix C. We then present the regularity assumptions and technical proofs for the main theoretical results in Appendices C and D, respectively. Finally, we include additional experimental results and further discussion of related work in Appendices G and F.

Appendix A More Details on Main Challenges and Contributions

Model-based clustering with variable selection, particularly using the sophisticated SRUW framework [38], offers a powerful approach to understanding data structure by not only grouping data but also by assigning distinct functional roles to variables. However, extending this rich framework to handle missing data, especially MNAR, and to perform variable selection efficiently in high-dimensional settings presents substantial theoretical and practical challenges. We address some difficulties as follows:

Integrating MNAR inside model-based clustering with variable selection: Unlike MAR or MCAR data, MNAR mechanisms are generally non-ignorable, meaning the missing data process must be explicitly modeled alongside the data distribution to avoid biased parameter estimates and incorrect clustering [33, 56]. Modeling this joint distribution significantly increases model complexity and the risk of misspecification. While various MNAR models exist (see [35, 62]), integration into complex variable selection clustering frameworks like SRUW, along with proving identifiability and estimator consistency, remains a frontier. Specifically, the MNARz mechanism (where missingness depends only on latent class membership) offers a tractable yet meaningful way to handle informative missingness, but its properties within a structured model like SRUW require careful elucidation.

Variable selection with SRUW and missing data in high-dimensional regime: The SRUW model’s strength lies in its detailed variable role specification. However, determining these roles from data is a combinatorial model selection problem. Original SRUW procedures often rely on stepwise algorithms that become computationally prohibitive as the number of variables D increases. Introducing missing data further complicates parameter estimation within each candidate model. Penalized likelihood methods, common in high-dimensional regression (e.g., LASSO [70, 75, 11]), offer a promising avenue for simultaneous parameter estimation and variable selection. However, adapting these to the GMM likelihood and then to the structured SRUW model with an integrated MNARz mechanism requires establishing theoretical guarantees such as parameter estimation consistency and ranking/selection consistency for the chosen penalized objective. Some notable, nontrivial hurdles are appropriate curvature properties like RSC for non-convex log-likelihood loss functions in high dimensions or the cone conditions to ensure sparse recovery.

Consistency of multi-stage procedures: Practical algorithms for variable selection tasks often involve multiple stages (e.g., initial ranking/variable screening, then detailed role assignment). Proving the consistency of such a multi-stage procedure requires advanced tools to handle.

Being aware of these difficulties, we attempt to make several contributions to address the aforementioned challenges, particularly focusing on providing a theoretical foundation for a two-stage variable selection procedure within the SRUW-MNARz framework:

  1. 1.

    Though the extension of the reinterpretation result from [62] is straightforward, this significantly simplifies the analysis and provides a way to handle identifiability of SRUW under MNAR missingness, which potentially facilitates algorithmic design.

  2. 2.

    Building on the MAR reinterpretation and existing results for SRUW data model identifiability and MNARz mechanism identifiability, we establish conditions under which the full SRUW-MNARz model parameters (including SRUW structure, GMM parameters, regression coefficients, and MNARz probabilities) are identifiable from the observed data. This is important since identifiability is a prerequisite for consistent parameter estimation.

  3. 3.

    To underpin the variable ranking stage, which uses a LASSO-like penalized GMM, we provide a rigorous analysis of the penalized GMM estimator based on [66, 44, 34]: First, we establish a non-asymptotic sup-norm bound on the gradient of the GMM negative log-likelihood at the true parameters, which is crucial for selecting the appropriate magnitude for regularization parameters. Subsequently, we formally state the RSC condition needed for the GMM loss function, and the error vector of the penalized GMM estimator lies in a specific cone, which is essential for sparse recovery. These settings provide the necessary high-dimensional statistical guarantees for the penalized GMM used in the ranking stage, ensuring that its estimated sparse means are reliable.

  4. 4.

    Leveraging the parameter consistency of the penalized GMM, we prove that the variable ranking score 𝒪K​(j) can consistently distinguish between truly relevant variables (for clustering means) and irrelevant variables. This involves establishing precise “signal strength” conditions; thereby, validating the first step of the two-stage variable selection procedure in [13].

  5. 5.

    For the main theorem on selection consistency, the process combines ranking consistency with an analysis of the BIC-based stepwise decisions for assigning roles. We show that under appropriate conditions on local BIC decision error rates and the choice of the stopping parameter c, the entire two-stage procedure consistently recovers the true variable partition. This provides the first (to our knowledge) formal consistency proof for such a two-stage variable selection in model-based clustering (SRUW) that incorporates an initial LASSO-like ranking, especially in the context of MNAR pattern. It bridges the theoretical gap between the computationally intensive and a more practical high-dimensional strategy.

  6. 6.

    While full optimality of c is not derived, we provide a condition for choosing c to achieve a target false positive rate in relevant variable selection, given the single-decision error probability. This moves beyond purely heuristic choices for c in [13], paving the way for more adaptive or theoretically grounded choices of c in practice. This offers a more principled approach to selecting this single hyperparameter in the algorithm, enhancing the reliability of the selection procedure.

Appendix B Preliminary of 𝕊​ℝ and 𝕊​ℝ​𝕌​𝕎 Models

This part summarizes definitions and theoretical results of the 𝕊​ℝ and 𝕊​ℝ​𝕌​𝕎 models in [37, 38]. We start with the following two definitions:

Definition 1 (𝕊​ℝ Model for Complete Data).

Let 𝐲n∈ℝD be the n-th observation, for n=1,…,N. A model ℳ𝕊​ℝ=(K,m,𝕊,ℝ) is defined by:

  1. 1.

    K: The number of mixture components (clusters).

  2. 2.

    m: The covariance structure of the GMM for variables in 𝕊.

  3. 3.

    𝕊⊆{1,…,D}: The non-empty set of relevant clustering variables.

  4. 4.

    ℝ⊆𝕊: The subset of relevant variables in 𝕊 used to explain the irrelevant variables 𝕊c={1,…,D}∖𝕊 via linear regression.

The parameters of the model are 𝛉=(𝛂,𝐚,𝛃,𝛀), where 𝛂 are the GMM parameters for variables in 𝕊. 𝐚 is the intercept vector, 𝛃 is the matrix of regression coefficients for explaining 𝐲𝕊c from 𝐲ℝ, and 𝛀 is the covariance matrix of the regression component. The complete data likelihood for one observation 𝐲n is:

f​(𝒚n;K,m,𝕊,ℝ,𝜽)=fclust​(𝒚n𝕊;K,m,𝜶)​freg​(𝒚n𝕊c;𝒂+𝒚nℝ​𝜷,𝛀)

where 𝐲n𝕊 denotes the sub-vector of 𝐲n corresponding to variables in 𝕊, and similarly for 𝐲n𝕊c and 𝐲nℝ. fclust is a K-component GMM density, and freg is a multivariate Gaussian density.

Importantly, 𝕊​ℝ model can be equivalently written as a single K-component GMM for the full data vector 𝐲n∈ℝD [37]:

f​(𝒚n;K,m,𝕊,ℝ,𝜽)=∑k=1Kπk​Φ​(𝒚n;𝝂k,𝚫k) (8)

where Φ​(⋅;𝛎,𝚫) is the multivariate Gaussian density with mean 𝛎 and covariance 𝚫. The parameters (𝛎k,𝚫k) are constructed from the original 𝕊​ℝ parameters (αk,𝐚,𝛃,𝛀) as follows: Let 𝚲 be a |𝕊|×|𝕊c| matrix derived from 𝛃 (which is |ℝ|×|𝕊c|). For a variable j∈𝕊 and l∈𝕊c, Λj​l=βj′​l if the j-th variable of 𝕊 is the j′-th variable of ℝ, and Λj​l=0 if the j-th variable of 𝕊 is in 𝕊∖ℝ. Then, for each component k:

  • •

    The mean vector 𝝂k∈ℝD has elements:

    νk​j={μk​jif variable ​j∈𝕊(𝒂+𝝁k𝕊​𝚲)jif variable ​j∈𝕊c

    where 𝝁k𝕊 is the mean of the k-th component for variables in 𝕊.

  • •

    The covariance matrix 𝚫k∈ℝD×D has blocks:

    Δk,j​l={Σk,j​lif ​j∈𝕊,l∈𝕊(𝚺k​𝚲)j​lif ​j∈𝕊,l∈𝕊c(𝚲T​𝚺k)j​lif ​j∈𝕊c,l∈𝕊(𝛀+𝚲T​𝚺k​𝚲)j​lif ​j∈𝕊c,l∈𝕊c
Definition 2 (𝕊​ℝ​𝕌​𝕎 Model for Complete Data).

Let 𝐲n∈ℝD be the n-th observation. An 𝕊​ℝ​𝕌​𝕎 model ℳ𝕊​ℝ​𝕌​𝕎=(K,m,r,l,𝕍) is defined by:

  1. 1.

    K: The number of mixture components.

  2. 2.

    m: The covariance structure of the GMM for variables in 𝕊.

  3. 3.

    r: The form of the covariance matrix 𝛀 for the regression of 𝕌 on ℝ.

  4. 4.

    l: The form of the covariance matrix 𝚪 for the independent variables 𝕎.

  5. 5.

    𝕍=(𝕊,ℝ,𝕌,𝕎): A partition of the D variables, where

    • •

      𝕊: Non-empty set of relevant clustering variables.

    • •

      ℝ⊆𝕊: Subset of 𝕊 regressing variables in 𝕌. (ℝ=∅ if 𝕌=∅).

    • •

      𝕌: Set of irrelevant variables explained by ℝ.

    • •

      𝕎: Set of irrelevant variables independent of all variables in 𝕊.

    • •

      𝕊∪𝕌∪𝕎={1,…,D} and these sets are disjoint.

The parameters are 𝛉=(𝛂,𝐚,𝛃,𝛀,𝛄,𝚪), where 𝛂 is the GMM parameters. The complete data density for one observation 𝐲n is:

f​(𝒚n;K,m,r,l,𝕍,𝜽)=fclust​(𝒚n𝕊;K,m,𝜶)​freg​(𝒚n𝕌;𝒂+𝒚nℝ​𝜷,𝛀)​findep​(𝒚n𝕎;𝜸,𝚪)

This can be written as a K-component GMM for the full data vector 𝐲n∈ℝD:

f​(𝒚n;K,m,r,l,𝕍,𝜽)=∑k=1Kπk​Φ​(𝒚n;𝝂~k,𝚫~k) (9)

where 𝛎~k and 𝚫~k are constructed from 𝛉. Let (𝛎k(𝕊,𝕌),𝚫k(𝕊,𝕌)) be the mean and covariance for the (𝕊,𝕌) part derived as in the 𝕊​ℝ model (where 𝕌 plays the role of 𝕊𝕊​ℝc). Then:

𝝂~k =((𝝂k(𝕊,𝕌))T,𝜸T)T
𝚫~k =(𝚫k(𝕊,𝕌)𝟎𝟎𝚪)

due to the independence of 𝕎 from 𝕊 and 𝕌 (given the cluster k, which is implicit in the construction of 𝛎k(𝕊,𝕌) using 𝛍k𝕊).

Similar to the global GMM representation of 𝕊​ℝ model under MAR, we can obtain such a global form of 𝕊​ℝ​𝕌​𝕎 model under the same missingness pattern, which is illustrated by the following claim.

Proposition 1.

Under MAR, 𝕊​ℝ​𝕌​𝕎 has a natural global GMM representation with appropriate parameters.

Proof of Proposition 1.

Let 𝒚n=(𝒚n𝕊,𝒚n𝕌,𝒚n𝕎) be the complete data vector for observation n∈{1,…,N}, partitioned according to the 𝕊​ℝ​𝕌​𝕎 model structure 𝕍=(𝕊,ℝ,𝕌,𝕎). The parameters are 𝜽=(𝜶,𝒂,𝜷,𝛀,𝜸,𝚪), where 𝜶=(π1,…,πk,𝝁1,…,𝝁K,𝚺1,…,𝚺K).

The complete-data density for a single observation 𝒚n, given its membership to cluster k (denoted by zi​k=1), is:

f​(𝒚n∣zi​k=1,𝕍,𝜽)=fclust​(𝒚n𝕊;𝝁k,𝚺k)⋅freg​(𝒚n𝕌;𝒂+𝒚nℝ​𝜷,𝛀)⋅findep​(𝒚n𝕎;𝜸,𝚪)

The marginal complete-data density for 𝒚n,by summing over latent cluster memberships is:

f​(𝒚n;𝕍,𝜽)=∑k=1Kπk​[fclust​(𝒚n𝕊;𝝁k,𝚺k)⋅freg​(𝒚n𝕌;𝒂+𝒚nℝ​𝜷,𝛀)⋅findep​(𝒚n𝕎;𝜸,𝚪)] (10)

Due to the independence of 𝒚n𝕎 from (𝒚n𝕊,𝒚n𝕌) given the cluster k and its parameters 𝜸,𝚪 are not k-specific, we can rewrite the term inside the sum as follows: suppose that fk​(𝒚n𝕊,𝒚n𝕌)=Φ​(𝒚n𝕊;𝝁k,𝚺k)⋅Φ​(𝒚n𝕌;𝒂+𝒚nℝ​𝜷,𝛀). This is the density of (𝒚n𝕊,𝒚n𝕌) given cluster k. Then, Equation 10 becomes:

f​(𝒚n;𝕍,𝜽)=(∑k=1Kπk​fk​(𝒚n𝕊,𝒚n𝕌))⋅findep​(𝒚n𝕎;𝜸,𝚪)

The term Φ​(𝒚n𝕎;𝜸,𝚪) factors out of the sum over k because 𝜸 and 𝚪 are not indexed by k.

Let 𝒚n=(𝒚no,𝒚nm) be the partition of 𝒚n into observed and missing parts for observation n. The missing parts are 𝒚nm=(𝒚n𝕊∩𝕄n,𝒚n𝕌∩𝕄n,𝒚n𝕎∩𝕄n), where 𝕄n denotes the set of missing variable indices for observation n. The observed-data density for observation n is obtained by integrating f​(𝒚n;𝕍,𝜽) over 𝒚nm:

f​(𝒚no;𝕍,𝜽) =∫f​(𝒚no,𝒚nm;𝕍,𝜽)​𝑑𝒚nm
=∫[(∑k=1Kπkfk(𝒚n𝕊∩𝕆n,𝒚n𝕊∩𝕄n,𝒚n𝕌∩𝕆n,𝒚n𝕌∩𝕄n))
×findep(𝒚n𝕎∩𝕆n,𝒚n𝕎∩𝕄n;𝜸,𝚪)]d𝒚n𝕊∩𝕄nd𝒚n𝕌∩𝕄nd𝒚n𝕎∩𝕄n

where 𝕆n denotes the set of observed variable indices for observation n. We can swap the sum and the integral:

f​(𝒚no;𝕍,𝜽) =∑k=1Kπk∫[fk(𝒚n𝕊∩𝕆n,𝒚n𝕊∩𝕄n,𝒚n𝕌∩𝕆n,𝒚n𝕌∩𝕄n)
×findep(𝒚n𝕎∩𝕆n,𝒚n𝕎∩𝕄n;𝜸,𝚪)]d𝒚n𝕊∩𝕄nd𝒚n𝕌∩𝕄nd𝒚n𝕎∩𝕄n

Since fk​(⋅) only involves variables in 𝕊 and 𝕌, and fclust​(⋅;𝜸,𝚪) only involves variables in 𝕎, and these sets are disjoint, the integration can be separated:

f​(𝒚no;𝕍,𝜽)=∑k=1Kπk [∫fk​(𝒚n𝕊∩𝕆n,𝒚n𝕊∩𝕄n,𝒚n𝕌∩𝕆n,𝒚n𝕌∩𝕄n)​𝑑𝒚n𝕊∩𝕄n​𝑑𝒚n𝕌∩𝕄n]
×[∫findep​(𝒚n𝕎∩𝕆n,𝒚n𝕎∩𝕄n;𝜸,𝚪)​𝑑𝒚n𝕎∩𝕄n]

Let fk​(𝒚n(𝕊,𝕌)∩𝕆n)=∫fk​(𝒚n𝕊∩𝕆n,𝒚n𝕊∩𝕄n,𝒚n𝕌∩𝕆n,𝒚n𝕌∩𝕄n)​𝑑𝒚n𝕊∩𝕄n​𝑑𝒚n𝕌∩𝕄n. This is the marginal density of the observed parts of (𝕊,𝕌) variables for component k. It is Gaussian, Φ​(𝒚n(𝕊,𝕌)∩𝕆n;𝝂k,o(𝕊,𝕌)​(n),𝚫k,o​o(𝕊,𝕌)​(n)), where these parameters are derived from the complete-data 𝕊​ℝ parameters for the (𝕊,𝕌) part. Let Φ​(𝒚n𝕎∩𝕆n;𝜸o(n),𝚪o​o(n))=∫Φ​(𝒚n𝕎∩𝕆n,𝒚n𝕎∩𝕄n;𝜸,𝚪)​𝑑𝒚n𝕎∩𝕄n. This is the marginal density of the observed parts of 𝕎 variables. Then,

f​(𝒚no;𝕍,𝜽)=(∑k=1Kπk​fclust​(𝒚n(𝕊,𝕌)∩𝕆n;𝝂k,o(𝕊,𝕌)​(n),𝚫k,o​o(𝕊,𝕌)​(n)))⋅findep​(𝒚n𝕎∩𝕆n;𝜸o(n),𝚪o​o(n)) (11)

From Equation 11, we can get

∏n=1N(∑k=1Kπk​fclust​(𝒚no;𝝂~k,o(n),𝚫~k,o​o(n)))

by allowing findep​(𝒚n𝕎∩𝕆n;𝜸o(n),𝚪o​o(n)) to be “absorbed” into the sum over k. So, we can write:

f​(𝒚no;𝕍,𝜽) =∑k=1K(πk​fclust​(𝒚n(𝕊,𝕌)∩𝕆n;𝝂k,o(𝕊,𝕌)​(n),𝚫k,o​o(𝕊,𝕌)​(n))⋅findep​(𝒚n𝕎∩𝕆n;𝜸o(n),𝚪o​o(n)))

Since the variables in (𝕊,𝕌) and 𝕎 are conditionally independent given k, the product of their marginal observed densities is the marginal observed density of their union. Let 𝒚no=(𝒚n(𝕊,𝕌)∩𝕆n,𝒚n𝕎∩𝕆n). Let 𝝂~k,o(n)=((𝝂k,o(𝕊,𝕌)​(n))T,(𝜸o(n))T)T. Let 𝚫~k,o​o(n)=(𝚫k,o​o(𝕊,𝕌)​(n)𝟎𝟎𝚪o​o(n)). Then, fclust​(𝒚n(𝕊,𝕌)∩𝕆n;𝝂k,o(𝕊,𝕌)​(n),𝚫k,o​o(𝕊,𝕌)​(n))⋅findep​(𝒚n𝕎∩𝕆n;𝜸o(n),𝚪o​o(n))=fclust​(𝒚no;𝝂~k,o(n),𝚫~k,o​o(n)).

Therefore, the observed-data density for a single observation n is:

f​(𝒚no;𝕍,𝜽)=∑k=1Kπk​fclust​(𝒚no;𝝂~k,o(n),𝚫~k,o​o(n)) (12)

And for N i.i.d. observations, the total observed-data likelihood is:

∏n=1Nf​(𝒚no;K,m,r,l,𝕍,𝜽)=∏n=1N(∑k=1Kπk​fclust​(𝒚no;𝝂~k,o(n),𝚫~k,o​o(n))) (13)

∎

When the parameters meet some distinct properties, 𝕊​ℝ​𝕌​𝕎 under complete-data can achieve identifiability, as shown in Theorem 1.

Fact 1 (Theorem 1 in [38]).

Let 𝚯(K,m,r,l,𝕍) be a subset of the parameter set Υ(K,m,r,l,𝕍) such that elements 𝛉=(𝛂,𝐚,𝛃,𝛀,𝛄,𝚪)

  • •

    contain distinct couples (μk,Σk) fulfilling ∀s⊆𝕊,∃(k,k′),1≤k<k′≤K;

    μk,s¯|s≠μk′,s¯|s​ or ​Σk,s¯|s≠Σk′,s¯|s​ or ​Σk,s¯​s¯|s≠Σk′,s¯​s¯|s (14)

    where s¯ denotes the complement in 𝕊 of any nonempty subset s of 𝕊

  • •

    if 𝕌≠∅,

    • –

      for all variables j of ℝ, there exists a variable u of 𝕌 such that the restriction 𝜷u​j of the regression coefficient matrix 𝜷 associated to j and u is not equal to zero.

    • –

      for all variables u of 𝕌, there exists a variable j of ℝ such that 𝜷u​j≠0.

  • •

    parameters Ω and τ exactly respect the forms r and l respectively. They are both diagonal matrices with at least two different eigenvalues if r=[L​B] and l=[L​B] and Ω has at least a non-zero entry outside the main diagonal if r=[L​C].

Let (K,m,r,l,𝕍) and (K∗,m∗,r∗,l∗,𝕍∗) be two models. If there exist 𝛉∈𝚯(K,m,r,l,𝕍) and 𝛉∗∈𝚯(K∗,m∗,r∗,l∗,𝕍∗) such that

f(⋅|K,m,r,l,𝕍,𝜽)=f(⋅|K∗,m∗,r∗,l∗,𝕍∗,𝜽∗)

then (K,m,r,l,𝕍)=(K∗,m∗,r∗,l∗,𝕍∗) and 𝛉=𝛉∗ (up to a permutation of mixture components).

Appendix C Regular Assumptions for Main Theoretical Results

For a fixed tuple of (K0,m0,r0,l0,𝕍0), we can write f​(𝒚;𝜽) instead of f​(𝒚;K0,m0,𝕊0,ℝ0,𝜽) or f​(𝒚;𝕍,𝜽) for short notation.

Assumption 1 (𝕊​ℝ).

There exists a unique (K0,m0,𝕊0,ℝ0) such that

h=f​(⋅;𝜽(K0,m0,𝕊0,ℝ0)∗)

for some parameter value 𝛉∗, and the pair (K0,m0) is assumed known. (In the following, we omit explicit notation of the dependence on (K0,m0) for brevity.)

Assumption 2.

The vectors 𝛉(𝕊,ℝ)∗ and 𝛉^(𝕊,ℝ) belong to a compact subset

𝚯(𝕊,ℝ)′⊂𝚯(𝕊,ℝ)

defined by

𝚯(𝕊,ℝ)′=𝒫K0−1×ℬ​(η,|𝕊|)K0×𝒟|𝕊|K0×ℬ​(ρ,|𝕊c|)×ℬ​(ρ,|ℝ|,|𝕊c|)×𝒟|𝕊c|

where the components are defined as follows:

  • •

    𝒫K−1={(π1,…,πK)∈[0,1]K:∑k=1Kπk=1} (set of proportions);

  • •

    ℬ​(η,r)={x∈ℝr:‖x‖≤η} with ‖x‖=∑i=1rxi2;

  • •

    ℬ​(ρ,q,r)={A∈ℳq×r​(ℝ):‖|A|‖≤ρ}, with norm

    ‖|A|‖=supy∈ℝr,‖y‖=1‖A​y‖;
  • •

    𝒟r is the set of r×r positive definite matrices with eigenvalues in [sm,sM] (with 0<sm<sM)

Assumption 3.

The optimal parameter 𝛉(𝕊0,ℝ0)∗ is an interior point of 𝚯(𝕊0,ℝ0)′.

The assumptions can be extended to the 𝕊​ℝ​𝕌​𝕎 framework. We define the optimal parameter as 𝜽𝕍∗, and the maximum likelihood estimator as 𝜽^𝕍. In this modified framework, the assumptions become:

Assumption 4 (𝕊​ℝ​𝕌​𝕎).

There exists a unique (K0,m0,r0,l0,𝕍0) such that

h=f​(⋅;𝜽(K0,m0,r0,l0,𝕍0)⋆)

for some 𝛉⋆, and the model (K0,m0,r0,l0,𝕍0) is assumed known.

Assumption 5.

Vectors 𝛉𝕍∗ and 𝛉^𝕍 belong to the compact subspace 𝚯𝕍′ of the compact set 𝚯𝕍, defined by:

𝚯𝕍′ =𝒫K0−1×ℬ​(η,|𝕊|)K0×𝒟|𝕊|K0×ℬ​(ρ,|𝕌|)×ℬ​(ρ,|ℝ|,|𝕌|)×𝒟|𝕌|×ℬ​(η,|𝕎|)×𝒟|𝕎|.
Assumption 6.

The optimal parameter 𝛉𝕍∗ is an interior point of 𝚯𝕍′.

Appendix D Proof of Main Results

Throughout the proofs, we denote the full parameter vector of the K-component Gaussian Mixture Model (GMM) by 𝜶=(𝝅,{𝝁k}k=1K,{𝚿k}k=1K), where 𝝅 are the mixing proportions, 𝝁k∈ℝD are the component means for D variables, and 𝚿k are the component precision matrices (inverse of covariance matrices 𝚺k). The true parameter vector is denoted by 𝜶∗. The negative single-observation log-likelihood is ℓ1​(𝒚n;𝜶)=−ln⁡fclust​(𝒚n;𝜶), and its sample average is ℓN​(𝜶). The score function components are Sj​(𝒚;𝜶∗)=−∂ℓ1​(𝒚;𝜶)∂αj|𝜶=𝜶∗. The Fisher Information Matrix (FIM) is 𝑰​(𝜶∗), and the empirical Hessian is 𝑯N​(𝜶). The set of true non-zero penalized parameters in 𝜶∗ is 𝕊0, with cardinality s0. The error vector is 𝚫=𝜶^−𝜶∗. For ranking, 𝒪K​(j) is the ranking score for variable j based on a grid of regularization parameters 𝒢λ. Standard notations like 𝔼​[⋅] for expectation, ∥⋅∥1,∥⋅∥2,∥⋅∥∞ for L1,L2,L∞ norms. Other abbreviations include KKT (Karush-Kuhn-Tucker), RSC (Restricted Strong Convexity), and SRUW (for the variable role framework).

D.1 Proof of Theorem 1

Proof Sketch. The proof involves two main steps. First, we establish that equality of observed-data likelihoods under MAR implies equality of the parameters of an equivalent full-data GMM representation for the 𝕊​ℝ​𝕌​𝕎 model. Second, we apply the identifiability arguments for the complete-data 𝕊​ℝ​𝕌​𝕎 model similar to [38], which itself relies on the identifiability of the 𝕊​ℝ model [37]

Proof.

Let ℳ1=(K,m,r,l,𝕍,𝜽) and ℳ2=(K∗,m∗,r∗,l∗,𝕍∗,𝜽∗) be two 𝕊​ℝ​𝕌​𝕎 models with ℳ1,ℳ2∈𝚯(K,m,r,l,𝕍). The complete-data density for an observation 𝒚n under ℳ1 and Equation 12 is given by:

f​(𝒚no;K,m,r,l,𝕍,𝜽,𝑴n)=∑k=1Kpk​fclust​(𝒚no;𝝂~k,o(n),𝚫~k,o​o(n)),

Analogously, we can derive the same expression for ℳ2. Consider the specific missingness pattern 𝑴0 where no data are missing for observation n. For this pattern, 𝒚no=𝒚n, 𝝂~k,o(n)=𝝂~k, and 𝚫~k,o​o(n)=𝚫~k. Under the premise of the theorem, if the equality of observed-data likelihoods for pattern 𝑴0, we get:

∑k=1Kπk​fclust​(𝒚n;𝝂~k,𝚫~k)=∑k′=1K∗pk′∗​fclust​(𝒚n;𝝂~k′∗,𝚫~k′∗),∀𝒚n∈ℝD.

From theorem 1 in [38], since 𝕊​ℝ​𝕌​𝕎 are identifiable under complete-data scenario this implies K=K∗, and up to a permutation of component labels:

πk=πk∗,𝝂~k=𝝂~k∗,and𝚫~k=𝚫~k∗for all ​k=1,…,K. (15)

For any 𝒚∈ℝD:

freg​(𝒚𝕌;r,𝒂+𝒚ℝ​𝜷,𝛀)​findep​(𝒚𝕎;l,𝜸,𝚪)=fclust​(𝒚𝕌∪𝕎;𝒂ˇ+𝒚ℝ​𝜷ˇ,𝛀ˇ)

where 𝕊c=𝕌∪𝕎, 𝒂ˇ=((𝒂)T,(𝜸)T)T, 𝜷ˇ=(𝜷𝟎) (with |𝟎|=|𝕎|×|ℝ|), and 𝛀ˇ=blockdiag​(𝛀,𝚪). Here, blockdiag​(𝛀,𝚪) is the block diagonal matrix created by aligning the input matrices 𝛀,𝚪. The 𝕊​ℝ​𝕌​𝕎 density can be expressed as an equivalent 𝕊​ℝ model density:

f​(𝒚;K,m,r,l,𝕍,𝜽) =fclust​(𝒚𝕊;K,m,𝜶)​fclust​(𝒚𝕊c;𝒂ˇ+𝒚ℝ​𝜷ˇ,𝛀ˇ)
=f𝕊​ℝ-eq​(𝒚;K,m,𝕊,ℝ,𝜽𝕊​ℝ-eq)

where 𝜽𝕊​ℝ-eq=(𝜶,𝒂ˇ,𝜷ˇ,𝛀ˇ). From Section D.1, the equivalent 𝕊​ℝ densities for the two models, where the second model is denoted with ∗) must be equal:

f𝕊​ℝ-eq​(𝒚;K,m,𝕊,ℝ,𝜽𝕊​ℝ-eq)=f𝕊​ℝ-eq​(𝒚;K∗,m∗,𝕊∗,ℝ∗,𝜽𝕊​ℝ-eq∗)

By the identifiability of the 𝕊​ℝ model [37] and three conditions in Theorem 1 , we deduce: K=K∗, m=m∗, 𝕊=𝕊∗, ℝ=ℝ∗, 𝜶=𝜶∗, 𝒂ˇ=𝒂ˇ∗, 𝜷ˇ=𝜷ˇ∗, and 𝛀ˇ=𝛀ˇ∗.

Since 𝕊=𝕊∗ and 𝕊c=𝕊∗c, supposed that, for contradiction, that 𝕌∗∩𝕎≠∅. Let j∈𝕌∗∩𝕎. Since j∈𝕎, for all q∈𝕌∗ the rows of (𝜷ˇ)q​j are 𝟎T. Additionally, since j∈𝕌∗, there exists some q∈ℝ∗=ℝ such that (𝜷ˇ∗)q​j≠𝟎T. However, we have established 𝜷ˇ=𝜷ˇ∗. This is a contradiction: (𝜷ˇ)q​j must be simultaneously 𝟎T and non-zero for all q∈s​R. Therefore, 𝕌∗∩𝕎=∅. By symmetry, 𝕌∩𝕎∗=∅. Given 𝕌∪𝕎=𝕌∗∪𝕎∗, the disjoint conditions imply 𝕌⊆𝕌∗ and 𝕌∗⊆𝕌, hence 𝕌=𝕌∗. Consequently, 𝕎=𝕎∗. This establishes 𝕍=𝕍∗. The parameter equalities 𝒂=𝒂∗,𝜷=𝜷∗,𝜸=𝜸∗,𝛀=𝛀∗,𝚪=𝚪∗ and r=r∗,l=l∗ then follow directly from comparing the components of 𝒂ˇ=𝒂ˇ∗, 𝜷ˇ=𝜷ˇ∗, 𝛀ˇ=𝛀ˇ∗, condition 3 in Theorem 1. Thus, the SRUW parameters under MAR are identifiable.

Now, we extend this to the full model which includes the MNARz mechanism parameters 𝝍={ψk}k=1K. The full set of parameters for the MNARz-SRUW model is (𝜽,𝝍). Suppose that the two models ℳA and ℳB produce the same observed-data likelihood for all observed data (𝒚no,𝒄n):

LMNARz-SRUW​(𝒚no,𝒄n;𝜽Adata,𝝍AMNARz)=LMNARz-SRUW​(𝒚no,𝒄n;𝜽Bdata,𝝍BMNARz)∀(𝒚no,𝒄n),

By Theorem 5, the equality of observed-data likelihoods under MNARz-SRUW:

LMNARz-SRUW​(𝒚no,𝒄n;𝜽Adata,𝝍AMNARz)=LMNARz-SRUW​(𝒚no,𝒄n;𝜽Bdata,𝝍BMNARz)

for all (𝒚no,𝒄n) implies the equality of the likelihoods for the augmented observed data 𝒚~no=(𝒚no,𝒄n) under the MAR interpretation where f​(𝒄n|zn​k=1;ψk) is treated as part of the component density:

∑k=1KA(πA)k​fk​(𝒚no|zn​k=1;(θASRUW)k)​f​(𝒄n|zn​k=1;(ψA)k)
=∑k′=1KB(πB)k′​fk′​(𝒚no|zn​k′=1;(θBSRUW)k′)​f​(𝒄n|zn​k′=1;(ψB)k′) (16)

for all (𝒚no,𝒄n). Let gk​(𝒚no,𝒄n;(θA)k,(ψA)k)=fk​(𝒚no|zn​k=1;(θASRUW)k)​f​(𝒄n|zn​k=1;(ψA)k). This means two finite mixture models on the augmented data (𝒚no,𝒄n) are identical. If ∑k=1KA(πA)k​gk​(⋅)=∑k′=1KB(πB)k′​gk′′​(⋅), and the family of component densities {gk} (and {gk′′}) is identifiable and satisfies certain linear independence conditions ([68], [2]), then KA=KB=K, and after a permutation of labels, (πA)k=(πB)k and gk​(⋅)=gk′​(⋅) for all k. Then, we need to ensure the family of component densities hk​(𝒚no,𝒄n;θk,ψk)=fk​(𝒚no|zn​k=1;θk)​f​(𝒄n|zn​k=1;ψk) meets such conditions. Hence, we further assume:

  1. 1.

    The parameter sets 𝜽SRUW and 𝝍MNARz are functionally independent

  2. 2.

    The support of 𝒚n and 𝒄n is rich enough such that equalities holding for all (𝒚n,𝒄n) imply functional equalities. This implies there is sufficient variability of 𝒚n and 𝒄n.

These assumptions hold because ψ parametrizes only f​(𝒄n|𝒛) and θ only f​(𝒚no|𝒛); hence there is no functional overlap. Applying those assumptions to Section D.1, we conclude that KA=KB=K, and up to a permutation of labels:

(πA)k=(πB)kfor all ​k, (17)

and

fk​(𝒚no|zn​k=1;(θASRUW)k)​f​(𝒄n|zn​k=1;(ψA)k)=fk​(𝒚no|zn​k=1;(θBSRUW)k)​f​(𝒄n|zn​k=1;(ψB)k) (18)

for all (𝒚no,𝒄n) and for each k=1,…,K.

Now we need to deduce that (θASRUW)k=(θBSRUW)k and (ψA)k=(ψB)k. Equation 18 can be written as:

A1​(𝒚no)​B1​(𝒄n)=A2​(𝒚no)​B2​(𝒄n)

where A1​(𝒚no)=fk​(𝒚no|zn​k=1;(θASRUW)k), B1​(𝒄n)=f​(𝒄n|zn​k=1;(ψA)k). By defined assumptions, this equality holding for all 𝒚no,𝒄n implies a relationship between these functions. Since 𝒚no and 𝒄n can vary independently 111Since, for each component k, the joint component density factorizes as fk​(𝒚no|zn​k=1;θk)​f​(𝒄n|zn​k=1;ψk) and both factors are positive on sets of positive measure, pick c0 in that common positive-measure set; integrating A1​(𝒚no)​B1​(c0)=A2​(𝒚no)​B2​(c0) over 𝒚no gives B1​(c0)=B2​(c0), hence A1=A2, and then B1=B2. and the parameters (θA)k,(ψA)k are functionally independent, if ∫A1​B1​𝑑𝒚o​𝑑𝒄=∫A2​B2​𝑑𝒚o​𝑑𝒄=1, and A1​B1=A2​B2 pointwise, and assuming A1,B1,A2,B2 are non-zero almost everywhere on their support: Pick a 𝒄n0 such that f​(𝒄n0|k;(ψA)k)≠0 and f​(𝒄n0|k;(ψB)k)≠0. Then for this fixed 𝒄n0:

fk​(𝒚no|zn​k=1;(θASRUW)k)⋅CA=fk​(𝒚no|zn​k=1;(θBSRUW)k)⋅CB

where CA=f​(𝒄n0|k;(ψA)k) and CB=f​(𝒄n0|k;(ψB)k). Since both fk​(𝒚no|…) are conditional densities , this implies CA=CB and thus

fk​(𝒚no|zn​k=1;(θASRUW)k)=fk​(𝒚no|zn​k=1;(θBSRUW)k)for all ​𝒚no.

By identifiability of SRUW under MAR, this implies that the structural SRUW parameters (mA,rA,lA,𝕍A) must equal (mB,rB,lB,𝕍B) and the parameters (θASRUW)k=(θBSRUW)k.

Since CA=CB means f​(𝒄n0|k;(ψA)k)=f​(𝒄n0|k;(ψB)k), and this must hold for all 𝒄n0 (by choosing different fixed 𝒄n0 or by varying 𝒚no first), this implies:

f​(𝒄n|k;(ψA)k)=f​(𝒄n|k;(ψB)k)for all ​𝒄n.

By idenitifiability of MNARz, this implies (ψA)k=(ψB)k.

Thus, we have KA=KB=K, and up to a common permutation of labels for k=1,…,K:

(πA)k=(πB)k
(mA,rA,lA,𝕍A)=(mB,rB,lB,𝕍B)
(θASRUW)k=(θBSRUW)k
(ψA)k=(ψB)k

This implies the full set of model specifications and parameters for ℳA and ℳB is identical. Therefore, the SRUW model under the specified MNARz mechanism is identifiable. ∎

D.2 Proof of Theorem 2

Proof Sketch.

Suppose that the observations are from a parametric model 𝒫𝚯={f​(⋅;𝜽):𝜽∈𝚯}, and assume that the distributions in 𝒫𝚯 are dominated by a common σ-finite measure υ with respect to which they have probability density functions f​(⋅;𝜽).

Now, we need to prove the consistency of the sample KL, which will lead to the consistency of the BIC criterion in both the 𝕊​ℝ and 𝕊​ℝ​𝕌​𝕎 models under the MAR mechanism. First, we prove the consistency of the sample KL in the 𝕊​ℝ model:

Proposition 2.

Under Assumption 1 and 2, for all (𝕊,ℝ),

1N​∑n=1Nln⁡(h​(𝒚no)f​(𝒚no;𝜽^(𝕊,ℝ)))​→n→∞𝑃​DKL​[h,f​(⋅;𝜽(𝕊,ℝ)∗)].
Proof.

Let 𝕆⊂{1,…,D} denote the indices of the observed components.

Recall that the full covariance matrix for the component k, 𝚫k, is built as:

Δk,j​l={Σk,j​l,j,l∈𝕊,(𝚺k​𝚲)j​l,j∈𝕊,l∈𝕊c,(𝚲T​𝚺k)j​l,j∈𝕊c,l∈𝕊,(𝛀+𝚲T​𝚺k​𝚲)j​l,j,l∈𝕊c.

For indices j,l∈𝕊: As before,

Δk,j​l=Σk,j​l.

Thus, eigenvalues of Δk,𝕊​𝕊 are in [sm,sM].

For j∈𝕊,l∈𝕊c:

Δk,j​l=(𝚺k​𝚲)j​l.

The norm of the product satisfies

∥𝚺k​𝚲∥≤∥𝚺k∥⋅∥𝚲∥.

Since 𝚺k∈𝒟|𝕊|,

∥𝚺k∥=λmax​(𝚺k)≤sM.

Given 𝜷∈ℬ​(ρ,|ℝ|,|𝕊c|),

∥𝚲∥≤ρ.

Thus,

∥𝚺k​𝚲∥≤sM​ρ.

This bounds the norm of the (𝕊×𝕊c) block of 𝚫k.

For j,l∈𝕊c:

Δk,j​l=(𝛀+𝚲T​𝚺k​𝚲)j​l.

We can bound the eigenvalues of this sum using the spectral norm:

∥𝛀+𝚲T​𝚺k​𝚲∥≤∥𝛀∥+∥𝚲T​𝚺k​𝚲∥.

Since 𝛀∈𝒟|𝕊c|,

∥𝛀∥=λmax​(𝛀)≤sM.

Next,

∥𝚲T​𝚺k​𝚲∥≤∥𝚲T∥⋅∥𝚺k∥⋅∥𝚲∥=∥𝚲∥2⋅∥𝚺k∥≤ρ2​sM.

Thus,

∥𝛀+𝚲T​𝚺k​𝚲∥≤sM+ρ2​sM=sM​(1+ρ2).

Given the above bounds, there exist refined constants s~m>0 and

s~M=max⁡{sM​(1+ρ2),sM​ρ,sM}=sM​(1+ρ2)

such that

s~m​I⪯𝚫k⪯s~M​I.

Note that s~M will depend on the sizes of the blocks, and it is finite being dependent on sM, ρ, and the structure of 𝛀, 𝜷, and 𝚺k.

Since 𝚫k,o​o is a principal submatrix of 𝚫k and 𝚫k is a symmetric matrix, by the interlacing theorem:

s~m≤λmin​(𝚫k,o​o)≤λmax​(𝚫k,o​o)≤s~M.

Thus,

s~m​I|𝕆|⪯𝚫k,o​o⪯s~M​I|𝕆|,

and

(𝚫k,o​o)−1⪯1s~m​I|𝕆|.

We are now moving to the main proof. By hypothesis, the true observed-data density is given by

h​(𝒚o)=∑k=1Kπk​Φ​(𝒚o;𝝂k,o∗,Δk,o​o∗)

and if we can show that

𝔼𝒚o​[|ln⁡h​(𝒚o)|]<∞.

then by the LLN,

1N​∑n=1Nln⁡(h​(𝒚no))​→n→∞𝑃​𝔼𝒚o​[ln⁡h​(𝒚o)].

Moreover, the observed-data likelihood under the model is

f​(𝒚o;𝜽)=∑k=1Kπk​Φ​(𝒚o;𝝂k,o,𝚫k,o​o).

and we later show in Proposition 3 that

1N​∑n=1Nln⁡(f​(𝒚no;𝜽^(𝕊,ℝ)))​→n→∞𝑃​𝔼𝒚o​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)].

To do so, we must verify that the class of functions

ℱ(𝕊,ℝ)={𝒚o↦ln⁡f​(𝒚o;𝜽):𝜽∈𝚯(𝕊,ℝ)′}

satisfies the conditions for uniform convergence under Assumption 2. In particular, by Assumption 2, the parameter space 𝚯(𝕊,ℝ)′ is compact, and for every 𝒚o∈ℝ|𝕆| the mapping

𝜽↦ln⁡f​(𝒚o;𝜽)

is continuous. Next, we verify that there is an h-integrable envelope function F∈ℱ(𝕊,ℝ). Recalling that for 𝒚o∈ℝ|𝕆|, the Gaussian component density is

Φ​(𝒚o;𝝂k,o,𝚫k,o​o)=(2​π)−|𝕆|2​|𝚫k,o​o|−12​exp⁡(−12​(𝒚o−𝝂k,o)T​𝚫k,o​o−1​(𝒚o−𝝂k,o)).

We can derive the bounds for this as follows.

For the upper bound, since λmin​(𝚫k,o​o)≥s~m,

|𝚫k,o​o|−12≤(s~m)−|𝕆|2.

With exp(−⋅)≤1 and 𝚫k,o​o is positive definite

Φ​(𝒚o;𝝂k,o,𝚫k,o​o)≤(2​π)−|𝕆|2​(s~m)−|𝕆|2.

Summing over k and using ∑k=1Kπk=1,

f​(𝒚o;𝜽)≤(2​π​s~m)−|𝕆|2.

Taking logarithms:

ln⁡f​(𝒚o;𝜽)≤−|𝕆|2​ln⁡(2​π​s~m).

For the lower bound, for each k,

ln⁡Φ​(𝒚o;𝝂k,o,𝚫k,o​o)=−|𝕆|2​ln⁡(2​π)−12​ln⁡|𝚫k,o​o|−12​(𝒚o−𝝂k,o)T​𝚫k,o​o−1​(𝒚o−𝝂k,o).

Since λmax​(𝚫k,o​o)≤s~M,

ln⁡|𝚫k,o​o|≤|𝕆|​ln⁡(s~M),

and 𝚫k,o​o−1⪯1s~m​I,

(𝒚o−𝝂k,o)T​𝚫k,o​o−1​(𝒚o−𝝂k,o)≤‖𝒚o−𝝂k,o‖2s~m.

Using ‖𝒚o−𝝂k,o‖2≤2​(‖𝒚o‖2+‖𝝂k,o‖2),

(𝒚o−𝝂k,o)T​𝚫k,o​o−1​(𝒚o−𝝂k,o)≤2​(‖𝒚o‖2+‖𝝂k,o‖2)s~m.

Given the construction of the mean vector 𝝂k for each component k:

νk​j={μk​j,if ​j∈𝕊,(a+𝝁k​Λ)j,if ​j∈𝕊c,

and knowing that a belongs to the closed ball ℬ​(ρ,1,|𝕊c|)-a set of 1×|𝕊c| matrices with norm bounded by ρ-we seek to derive a uniform bound for 𝝂k. For indices j∈𝕊, νk​j=μk​j, and since the parameter space 𝚯V′ is compact, the cluster means μk​j are bounded by some constant η>0. For indices j∈𝕊c, we have νk​j=(a+𝝁k​Λ)j=aj+(𝝁k​Λ)j. Using the elementary inequality (u+v)2≤2​(u2+v2), it follows that

∥νk​j∥2≤2​(∥aj∥2+∥(𝝁k​Λ)j∥2).

Given that a lies in ℬ​(ρ,1,|𝕊c|), ∥aj∥2≤ρ2 for each j∈𝕊c. In addition, the term (𝝁k​Λ)j can be bounded by ‖𝝁k‖2​‖Λ⋅j‖2, where Λ⋅j is the j-th column of Λ. With 𝝁k belonging to a compact set ℬ​(η,|𝕊|), ∥𝝁k∥2≤η2, and because Λ is derived from bounded parameters (including 𝜷 with norm being bounded by ρ2), each column Λ⋅j is bounded in norm by ρ2. Consequently, |(𝝁k​Λ)j|≤ρ2+η2​ρ2=ρ2​(1+η2) for each j∈𝕊c. Combining these results, for any j∈𝕊c,

∥νk​j∥2≤2​ρ2​(1+η2).

Thus, all entries of 𝝂k, irrespective of whether they correspond to indices in 𝕊 or 𝕊c, are bounded by a constant that depends on η, ρ. Hence,

‖𝝂k‖2 =∑j‖𝝂k,j‖2
≤∑j2​ρ2​(1+η2)
=2​D​ρ2​(1+η2)

uniformly for all k.

Consider again the Gaussian density for the observed variables

ln⁡Φ​(𝒚o;𝝂k,o,𝚫k,o​o)=−|𝕆|2​ln⁡(2​π)−12​ln⁡|𝚫k,o​o|−12​(𝒚o−𝝂k,o)T​𝚫k,o​o−1​(𝒚o−𝝂k,o).

To control the quadratic term in the exponent, we note that

‖𝒚o−𝝂k,o‖2≤2​(‖𝒚o‖2+‖𝝂k,o‖2).

Given that ‖𝝂k‖=∑j∈𝕆∥νk​j∥2≤2​|𝕆|​ρ2​(1+η2) is uniformly bounded, it follows that

‖𝒚o−𝝂k,o‖2≤2​(‖𝒚o‖2+2​|𝕆|​ρ2​(1+η2)).

This bound ensures that the quadratic form (𝒚o−𝝂k,o)T​𝚫k,o​o−1​(𝒚o−𝝂k,o) remains finite. The uniform boundedness of 𝝂k thus allows us to assert the existence of an integrable envelope function F​(𝒚o) that uniformly bounds |ln⁡f​(𝒚o;𝜽)| for all 𝜽 within the compact parameter space 𝚯V′.

Thus,

ln⁡Φ​(𝒚o;𝝂k,o,𝚫k,o​o)≥−|𝕆|2​ln⁡(2​π)−|𝕆|2​ln⁡(s~M)−12⋅2​(‖𝒚o‖2+2​|𝕆|​ρ2​(1+η2))s~m.

Simplifying,

ln⁡Φ​(𝒚o;𝝂k,o,𝚫k,o​o)≥−|𝕆|2​ln⁡(2​π​s~M)−∥𝒚o∥2+2|𝕆|ρ2(1+η2))s~m.

Using Jensen’s inequality over the mixture:

ln⁡f​(𝒚o;𝜽)≥∑k=1Kπk​(−|𝕆|2​ln⁡(2​π​s~M)−∥𝒚o∥2+2|𝕆|ρ2(1+η2))s~m).

Since ∑k=1Kπk=1,

ln⁡f​(𝒚o;𝜽)≥−|𝕆|2​ln⁡(2​π​s~M)−∥𝒚o∥2+2|𝕆|ρ2(1+η2))s~m.

Combining the refined upper and lower bounds:

−|𝕆|2​ln⁡(2​π​s~M)−∥𝒚o∥2+2|𝕆|ρ2(1+η2))s~m≤ln⁡f​(𝒚o;𝜽(𝕊,ℝ))≤−|𝕆|2​ln⁡(2​π​s~m).

As a final note, these bounds also rely on the eigenvalue constraints on 𝚺k, 𝛀, and the norm constraint on 𝜷.

Now we prove that the envelope function F, which is related to the upper and lower bounds above, is h-integrable. The true density h​(𝒚) corresponds to a Gaussian mixture model f​(𝒚;𝜽(𝕊0,ℝ0)∗) with parameters in a compact set. The observed-data density for 𝒚o, obtained by marginalizing over the missing components, is given by

h​(𝒚o)=∑k=1Kπk​Φ​(𝒚o;𝝂k,o∗,Δk,o​o∗).

To verify that the envelope function F is h-integrable, we need to show

∫‖𝒚o‖2​h​(𝒚o)​𝑑𝒚o<∞.

We proceed by examining the second moment of the observed components. First, observe that

∫‖𝒚o‖2​h​(𝒚o)​𝑑𝒚o =∑k=1Kπk​∫‖𝒚o‖2​Φ​(𝒚o;𝝂k,o∗,Δk,o​o∗)​𝑑𝒚o.
≤2​‖𝝂k,o∗‖2+2​tr⁡(Δk,o​o∗).

by using Lemma. By the compactness of the parameter space 𝚯𝕍 and since ∑k=1Kπk=1 and the bounds derived earlier, we have

∫‖𝒚o‖2​h​(𝒚o)​𝑑𝒚o≤|𝕆|​ρ2​(1+η2)+2​sM​|𝕊0,o|<∞

Therefore, F is h-integrable, i.e.,

∫F​(𝒚o)​h​(𝒚o)​𝑑𝒚o<∞,

since ∫‖𝒚o‖2​h​(𝒚o)​𝑑𝒚o<∞.

Finally, because |ln⁡(h​(𝒚o))|≤F​(𝒚o), it implies

𝔼​[|ln⁡(h​(𝒚o))|]=∫|ln⁡(h​(𝒚o))|​h​(𝒚o)​𝑑𝒚o≤∫F​(𝒚o)​h​(𝒚o)​𝑑𝒚o<∞.

Hence, the envelope function F is h-integrable and 𝔼​[|ln⁡(h​(𝒚o))|]<∞. Thus, we can apply the law of large numbers to obtain the consistency of the sample KL divergence. ∎

Proposition 3.

Assume that

  1. 1.

    (𝒚1o,…,𝒚no) are i.i.d. observed vectors with unknown density h.

  2. 2.

    𝚯 is a compact metric space.

  3. 3.

    𝜽∈𝚯↦ln⁡[f​(𝒚o;𝜽)] is continuous for every 𝒚o∈ℝ|𝕆|.

  4. 4.

    F is an envelope function of ℱ:={ln⁡[f​(⋅;𝜽)];𝜽∈𝚯} which is h-integrable.

  5. 5.

    𝜽∗=argmax𝜽∈𝚯DKL​[h,f​(⋅;𝜽)]

  6. 6.

    𝜽^=argmax𝜽∈𝚯​∑n=1Nln⁡f​(𝒚n;𝜽).

Then, as n→∞,

1N​∑n=1Nln⁡f​(𝒚no;𝜽^(𝕊,ℝ))​→n→∞𝑃​𝔼𝒚o​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)].
Proof.

We consider the following inequality:

|𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)]−1N​∑n=1Nln⁡f​(𝒚no;𝜽^(𝕊,ℝ))|
≤|𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)]−𝔼​[ln⁡f​(𝒚o;𝜽^(𝕊,ℝ))]|+sup𝜽∈𝚯(𝕊,ℝ)|𝔼​[ln⁡f​(𝒚o;𝜽)]−1N​∑n=1Nln⁡f​(𝒚no;𝜽)|.

By the definition of 𝜽(𝕊,ℝ)∗, we have

𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)]−𝔼​[ln⁡f​(𝒚o;𝜽^(𝕊,ℝ))]≥0.

Note that,

𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)]−𝔼​[ln⁡f​(𝒚o;𝜽^(𝕊,ℝ))] =𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)]−1N​∑n=1Nln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)⏟I
+1N​∑n=1Nln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)−1N​∑n=1Nln⁡f​(𝒚o;𝜽^(𝕊,ℝ))⏟I​I
−1N∑n=1Nlnf(𝒚o;𝜽^(𝕊,ℝ))+𝔼[lnf(𝒚o;𝜽^(𝕊,ℝ))⏟I​I​I]

For term I​I, sine 𝜽^ maximizes the empirical log-likelihood,

1N∑n=1Nln(f(Xn|𝜽^(𝕊,ℝ))≥1N∑n=1Nln(f(Xn|𝜽(𝕊,ℝ)∗),

which implies that I​I≤0. Therefore,

𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)]−𝔼​[ln⁡f​(𝒚o;𝜽^(𝕊,ℝ))]≤|I|+|I​I​I|.

Moreover, since both |I|​ and ​|I​I​I| are bounded by the uniform deviation:

|I|≤sup𝜽∈𝚯(𝕊,ℝ)|𝔼​[ln⁡f​(𝒚o;𝜽)]−1N​∑n=1Nln⁡f​(𝒚no;𝜽)|
and
|I​I|≤sup𝜽∈𝚯(𝕊,ℝ)|𝔼​[ln⁡f​(𝒚o;𝜽)]−1N​∑n=1Nln⁡f​(𝒚no;𝜽)|

Hence,

𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)]−𝔼​[ln⁡f​(𝒚o;𝜽^(𝕊,ℝ))] ≤|I|+|I​I​I|
≤2​sup𝜽∈𝚯(𝕊,ℝ)|𝔼​[ln⁡f​(𝒚o;𝜽)]−1N​∑n=1Nln⁡f​(𝒚no;𝜽)|.

Hence, the left-hand side is bounded by three times the uniform deviation:

|𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)]−1N​∑n=1Nln⁡f​(𝒚no;𝜽^(𝕊,ℝ))|≤3​sup𝜽∈𝚯(𝕊,ℝ)|𝔼​[ln⁡f​(𝒚o;𝜽)]−1N​∑n=1Nln⁡f​(𝒚no;𝜽)|.

By the argument in the Proposition 2-namely, using the compactness of 𝚯(𝕊,ℝ), the continuity of 𝜽↦ln⁡f​(⋅;𝜽), and the existence of an h-integrable envelope F-the class ℱ(𝕊,ℝ) is P-Glivenko-Cantelli by applying Example 19.8 in van der Vaart (1998) to conclude on the finiteness of bracketing numbers of ℱ under the assumptions. In particular,

sup𝜽∈𝚯(𝕊,ℝ)|𝔼​[ln⁡f​(𝒚o;𝜽)]−1N​∑n=1Nln⁡f​(𝒚no;𝜽)|​→n→∞𝑃​0.

Therefore,

1N​∑n=1Nln⁡f​(𝒚no;𝜽^(𝕊,ℝ))​→n→∞𝑃​𝔼​[ln⁡f​(𝒚o;𝜽(𝕊,ℝ)∗)].

which concludes the proof of Proposition 2. ∎

Now, we prove the consistency of the sample KL in the 𝕊​ℝ​𝕌​𝕎 model:

Proposition 4.

Under Assumption 4 and 5, for all 𝕍,

1N​∑n=1Nln⁡(h​(𝒚no)f​(𝒚no;𝜽^𝕍))​→n→∞𝑃​DKL​[h,f​(⋅;𝜽𝕍∗)].
Proof.

We carry the proof in the same manner as in the 𝕊​ℝ model. Let 𝕍=(𝕊,ℝ,𝕌,𝕎) and 𝕆⊂{1,…,D} denotes the observed variables. We still want to apply the Proposition 3 with the new family

ℱ(𝕍)≔{ln[f(⋅|𝜽)];𝜽∈𝚯𝕍′}

similarly to the proof of Proposition 2 to achieve

1N​∑n=1Nln⁡(f​(𝒚no;𝜽^𝕍))​→n→∞𝑃​𝔼𝒚o​[ln⁡f​(𝒚o;𝜽𝕍∗)].

To do so, we must initially verify that the class of functions

ℱ(𝕍)={𝒚o↦ln⁡f​(𝒚o;𝜽):𝜽∈𝚯𝕍′}

satisfies the conditions for uniform convergence under Assumption 5. In particular, by Assumption 5, the parameter space 𝚯𝕍′ is compact, and for every 𝒚o∈ℝ|𝕆| the mapping

𝜽↦ln⁡f​(𝒚o;𝜽)

is continuous. Next, we verify that there is an h-integrable envelope function F∈ℱ(𝕍). Recalling that for 𝒚o∈ℝ|𝕆| and a given component k, the current Gaussian component density is

Φ​(𝒚o;𝝂~k,o,𝚫~k,o​o)=(2​π)−|𝕆|2​|𝚫~k,o​o|−12​exp⁡(−12​(𝒚o−𝝂~k,o)T​𝚫~k,o​o−1​(𝒚o−𝝂~k,o)).

We will bound the density function as usual. From Proposition 2, we know that there exists s~m>0 and s~M′=max⁡{sM​(1+ρ2),sM​ρ,sM}=sM​(1+ρ2) such that

s~m​I⪯𝚫k⪯s~M​I.

Thus,

s~m​I|𝕆|⪯𝚫k,o​o⪯s~M​I|𝕆|,

and

(𝚫k,o​o)−1⪯1s~m​I|𝕆|.

Given that 𝚪∈𝒟|𝕎| and is independent of any relevant variables, we have that the principal sub-matrix 𝚪o​o is bounded by sM

λmax​(𝚪o​o)≤sM.

Since 𝚫~k,o​o is block-diagonal, its eigenvalues are the union of the eigenvalues of 𝚫k,o​o and 𝚪o​o. Hence,

λmax​(𝚫~k,o​o)≤max⁡{sM​(1+ρ2),sM}=sM​(1+ρ2).

We define the upper bound by

s~M:=sM​(1+ρ2).

The lower bound of the structured block, together with the lower bound on the independent block, implies that the covariance 𝚫~k,o​o satisfies

s~m​I⪯𝚫~k,o​o⪯s~M​I.

Because 𝚫~k,o​o−1⪯1s~m​I, it follows that

(𝒚o−𝝂~k,o)T​𝚫~k,o​o−1​(𝒚o−𝝂~k,o)≤‖𝒚o−𝝂~k,o‖2s~m.

Using the elementary inequality

‖𝒚o−𝝂~k,o‖2≤2​(‖𝒚o‖2+‖𝝂~k,o‖2),

we obtain

(𝒚o−𝝂~k,o)T​𝚫~k,o​o−1​(𝒚o−𝝂~k,o)≤2​(‖𝒚o‖2+‖𝝂~k,o‖2)s~m.

Now we process the bound for 𝝂~k,o. Recalling that if we denote D𝕊​𝕌 denote the number of coordinates in the structured part (i.e. 𝕊∪𝕌); then

‖𝝂k,o‖2=∑j∈o∥νk​j∥2≤D𝕊​𝕌,o​ρ2​(1+η2)

is uniformly bounded and ∥γo∥2≤η2 since γ belongs to the compact set ℬ​(η,|𝕎|), it follows that

‖𝝂~k,o‖2≤D𝕊​𝕌,o​ρ2​(1+η)2+η2,

We deduce that

ln⁡Φ​(𝒚o;𝝂~k,o,𝚫~k,o​o)≥−|𝕆|2​ln⁡(2​π​s~M)−‖𝒚o‖2+D𝕊​𝕌,o​ρ2​(1+η2)+η2s~m.

Using Jensen’s inequality over the mixture and ∑k=1Kπk=1, we have

ln⁡f​(𝒚o;𝜽) ≥∑k=1Kπk​(−|𝕆|2​ln⁡(2​π​s~M)−‖𝒚o‖2+D𝕊​𝕌,o​ρ2​(1+η2)+η2s~m).
≥−|𝕆|2​ln⁡(2​π​s~M)−‖𝒚o‖2+D𝕊​𝕌,o​ρ2​(1+η2)+η2s~m

Therefore, each family’s member in ℱ(𝕍) is bounded by

−|𝕆|2​ln⁡(2​π​s~M)−‖𝒚o‖2+D𝕊​𝕌,o​ρ2​(1+η2)+η2s~m≤ln⁡f​(𝒚o;𝜽(𝕍))≤−|𝕆|2​ln⁡(2​π​s~m).

Now we prove that F is h-integrable. The true observed-data density under the 𝕊​ℝ​𝕌​𝕎 model is given by

h​(𝒚o)=∑k=1Kπk​Φ​(𝒚o;𝝂~k,o∗,Δ~k,o​o∗),

where the “∗” indicates that the parameters are the true ones and the density is that of a Gaussian mixture with parameters in a compact set. In our derivation we have shown that, for any 𝜽 in the family ℱ(𝕍), the log-density satisfies

−|𝕆|2​ln⁡(2​π​s~M)−‖𝒚o‖2+D𝕊​𝕌,o​ρ2​(1+η2)+η2s~m≤ln⁡f​(𝒚o;𝜽)≤−|𝕆|2​ln⁡(2​π​s~m).

In other words, every function ln⁡f​(𝒚o;𝜽) is bounded in absolute value by

F​(𝒚o)=|𝕆|2​ln⁡(2​π​s~M)+‖𝒚o‖2+D𝕊​𝕌,o​ρ2​(1+η2)+η2s~m.

Hence, for all 𝜽 in the compact parameter space,

|ln⁡f​(𝒚o;𝜽)|≤F​(𝒚o).

To verify that F is h-integrable, we must show that

∫F​(𝒚o)​h​(𝒚o)​𝑑𝒚o<∞.

Since the envelope function F is of the form

F​(𝒚o)=C0+1s~m​‖𝒚o‖2,

with

C0=|𝕆|2​ln⁡(2​π​s~M)+D𝕊​𝕌,o​ρ2​(1+η2)+η2s~m,

it is clear that

∫F​(𝒚o)​h​(𝒚o)​𝑑𝒚o≤C0+1s~m​∫‖𝒚o‖2​h​(𝒚o)​𝑑𝒚o.

Because h​(𝒚o) is a Gaussian mixture with parameters in a compact set, standard properties of Gaussian mixtures guarantee that the second moment is finite, that is,

∫‖𝒚o‖2​h​(𝒚o)​𝑑𝒚o<∞.

Thus,

∫F​(𝒚o)​h​(𝒚o)​𝑑𝒚o<∞.

Finally, since for every 𝒚o we have

|ln⁡h​(𝒚o)|≤F​(𝒚o),

It follows that

𝔼​[|ln⁡h​(𝒚o)|]=∫|ln⁡h​(𝒚o)|​h​(𝒚o)​𝑑𝒚o≤∫F​(𝒚o)​h​(𝒚o)​𝑑𝒚o<∞.

Therefore, the envelope function F is h-integrable and 𝔼​[|ln⁡(h​(𝒚o))|]<∞, and consequently, we can apply the law of large numbers and uniform convergence to conclude the proof. ∎

Proof of Theorem 2.

Define the BIC score of a variable partition 𝕍=(𝕊,ℝ,𝕌,𝕎) by

BIC​(𝕍)=2​ℓN​(𝜽^𝕍)−Ξ(𝕍)​log⁡N,

where Ξ(𝕍) is the number of free parameters. Let 𝕍0 be the true partition and set Δ​BIC​(𝕍)=BIC​(𝕍0)−BIC​(𝕍).

Let 𝕍1={𝕍≠𝕍0:DKL​[h,f​(⋅;𝜽𝕍⋆)]>0}. The Identifiability Theorem 1 implies every 𝕍≠𝕍0 belongs to 𝕍1. For 𝕍∈𝕍1,

Δ​BIC​(𝕍)=2​N​[1N​∑n=1Nln⁡(f​(𝒚no;𝜽^𝕍0)h​(𝒚no))−1N​∑n=1Nln⁡(f​(𝒚no;𝜽^𝕍)h​(𝒚no))]+[Ξ(𝕍)−Ξ(𝕍0)]​log⁡N (19)

To prove the theorem, we will prove that:

∀𝕍∈𝕍1,ℙ​[Δ​BIC​(𝕍)<0]​→N→∞​0

Denoting 𝕄N​(𝕍)=1N​∑n=1Nln⁡(f​(𝒚no;𝜽^𝕍)h​(𝒚no)),M​(𝕍)=−DKL​[h,f​(⋅;𝜽𝕍∗)], from Equation 19, we get:

ℙ​(Δ​BIC​(𝕍)<0)
= ℙ​(2​N​(𝕄N​(𝕍0)−𝕄N​(𝕍))+[Ξ(𝕍)−Ξ(𝕍0)]​log⁡N<0)
= ℙ​(𝕄N​(𝕍0)−M​(𝕍0)+M​(𝕍0)−M​(𝕍)+M​(𝕍)−𝕄N​(𝕍)+[Ξ(𝕍)−Ξ(𝕍0)]​log⁡N2​N<0)
≤ ℙ​(𝕄N​(𝕍0)−M​(𝕍0)>ϵ)+ℙ​(M​(𝕍)−𝕄N​(𝕍)>ϵ)+ℙ​(M​(𝕍0)−M​(𝕍)+[Ξ(𝕍)−Ξ(𝕍0)]​log⁡N2​N<2​ϵ)

From Proposition 4, we know that

1N​∑n=1Nln⁡(h​(𝒚no)f​(𝒚no;𝜽^𝕍))​→n→∞𝑃​DKL​[h,f​(⋅;𝜽𝕍∗)].

This leads to 𝕄N​(𝕍)​→n→∞𝑃​M​(𝕍). Similar to the proof for the 𝕊​ℝ​𝕌​𝕎 model in Maugis [38], we prove the theorem. ∎

D.3 Proof of Theorem 3

Throughout this proof, we analyze the estimator defined by minimizing the negative log-likelihood penalized. Accordingly, although we denote by

ℓ​(𝒚;𝜶)=ln⁡f​(𝒚;𝜶)

the usual log-likelihood of a single observation, we consistently work with its negative −ℓ as the optimization objective. To streamline notation, we therefore use the term “score” to mean the gradient of the negative log-likelihood:

Sj​(𝒚n;𝜶):=−∂ℓ​(𝒚n;𝜶)∂αj.

This convention differs by a minus sign from the standard statistical score, but it ensures that all gradients, Hessians, and Fisher information matrices below are taken with respect to the minimized objective −ℓN.

We start with presenting some useful lemmas before proving our main theorem. These lemmas establish fundamental properties of the penalized GMM estimators and the variable ranking procedure.

Lemma 1 (Score Function Components).

Let fclust​(𝐲n;𝛂)=∑m=1Kπm​Φ​(𝐲n;𝛍m,𝚺m) be the p.d.f of a K-component GMM for an observation 𝐲n, where 𝛂=(𝛑,{𝛍k}k=1K,{𝚺k}k=1K) represents the GMM parameters. Let 𝚿k=𝚺k−1 be the precision matrix for component k. The log-likelihood for observation 𝐲n is ℓ​(𝐲n;𝛂)=ln⁡fclust​(𝐲n;𝛂). The responsibility for component k given observation 𝐲n and parameters 𝛂 is tk​(𝐲n;𝛂)=πk​Φ​(𝐲n;𝛍k,𝚺k)fclust​(𝐲n;𝛂). The score components, Sj​(𝐲n;𝛂∗), where αj is a parameter in 𝛂 and 𝛂∗ are the true parameter values, are given as follows:

  1. 1.

    Mixing proportions πk: When parameterizing πK=1−∑j=1K−1πj, the score component for πk, k∈{1,…,K−1}, is:

    Sπk​(𝒚n;𝜶∗)=tk​(𝒚n;𝜶∗)πK∗−tk​(𝒚n;𝜶∗)πk∗
  2. 2.

    Mean parameters (𝝁k)d: Let (𝝁k)d be the d-th component of the mean vector 𝝁k. The score component is:

    S(𝝁k)d​(𝒚n;𝜶∗)=−tk​(𝒚n;𝜶∗)​(𝚿k∗​(𝒚n−𝝁k∗))d
  3. 3.

    Precision matrix parameters (Ψk)r​s: Let (𝚿k)r​s be the element (r,s) of the symmetric precision matrix 𝚿k. The score component is:

    S(𝚿k)r​s​(𝒚n;𝜶∗)=−tk​(𝒚n;𝜶∗)​(2−δr​s)2​((𝚺k∗)r​s−(𝒚n−𝝁k∗)r​(𝒚n−𝝁k∗)s)

    where δr​s is the Kronecker delta.

Proof of Lemma 1.

Let Φm​(𝒚n;𝜶)=Φ​(𝒚n;𝝁m,𝚺m). The log-likelihood for a single observation 𝒚n is ℓ​(𝒚n;𝜶)=ln⁡(∑m=1Kπm​Φm​(𝒚n;𝜶)). The score component for a generic parameter αj is Sj​(𝒚n;𝜶∗)=−∂ℓ​(𝒚n;𝜶)∂αj|𝜶=𝜶∗.

Score components for mixing proportions πk. We parameterize πK=1−∑j=1K−1πj. For k∈{1,…,K−1}:

∂ℓ​(𝒚n;𝜶)∂πk =1fclust​(𝒚n;𝜶)​∂∂πk​(∑m=1K−1πm​Φm​(𝒚n;𝜶)+(1−∑j=1K−1πj)​ΦK​(𝒚n;𝜶))
=Φk​(𝒚n;𝜶)−ΦK​(𝒚n;𝜶)fclust​(𝒚n;𝜶)

Since tm​(𝒚n;𝜶)=πm​Φm​(𝒚n;𝜶)fclust​(𝒚n;𝜶), we have Φm​(𝒚n;𝜶)fclust​(𝒚n;𝜶)=tm​(𝒚n;𝜶)πm. Thus,

∂ℓ​(𝒚n;𝜶)∂πk=tk​(𝒚n;𝜶)πk−tk​(𝒚n;𝜶)πK

The score component at 𝜶∗ is:

Sπk​(𝒚n;𝜶∗)=−(tk​(𝒚n;𝜶∗)πk∗−tk​(𝒚n;𝜶∗)πK∗)=tk​(𝒚n;𝜶∗)πK∗−tk​(𝒚n;𝜶∗)πk∗.

Score components for mean parameters (μk)d. Let (𝝁k)d be the d-th component of 𝝁k.

∂ℓ​(𝒚n;𝜶)∂(𝝁k)d =1fclust​(𝒚n;𝜶)​∂∂(𝝁k)d​(∑m=1Kπm​Φm​(𝒚n;𝜶))
=πkfclust​(𝒚n;𝜶)​∂Φk​(𝒚n;𝜶)∂(𝝁k)d

Since ∂Φk∂(𝝁k)d=Φk​∂ln⁡Φk∂(𝝁k)d, and for Φk, ln⁡Φk​(𝒚n;𝜶)=Ck−12​(𝒚n−𝝁k)⊤​𝚿k​(𝒚n−𝝁k), we have:

∂ln⁡Φk​(𝒚n;𝜶)∂(𝝁k)d=(𝚿k​(𝒚n−𝝁k))d

Substituting this back:

∂ℓ​(𝒚n;𝜶)∂(𝝁k)d=πk​Φk​(𝒚n;𝜶)fclust​(𝒚n;𝜶)​(𝚿k​(𝒚n−𝝁k))d=tk​(𝒚n;𝜶)​(𝚿k​(𝒚n−𝝁k))d

The score component at 𝜶∗ is:

S(𝝁k)d​(𝒚n;𝜶∗)=−tk​(𝒚n;𝜶∗)​(𝚿k∗​(𝒚n−𝝁k∗))d.

Score components for precision matrix parameters (𝚿k)r​s. Let (𝚿k)r​s be an element of the symmetric precision matrix 𝚿k.

∂ℓ​(𝒚n;𝜶)∂(𝚿k)r​s=πkfclust​(𝒚n;𝜶)​∂Φk​(𝒚n;𝜶)∂(𝚿k)r​s=tk​(𝒚n;𝜶)​1πk​πk​∂ln⁡Φk​(𝒚n;𝜶)∂(𝚿k)r​s

The log-density of a single Gaussian component is ln⁡Φk​(𝒚n;𝜶)=Ck′+12​ln​det(𝚿k)−12​(𝒚n−𝝁k)⊤​𝚿k​(𝒚n−𝝁k). For a symmetric matrix 𝚿k, the derivative with respect to an element (𝚿k)r​s (using the “symmetric” derivative convention where ∂∂Xr​s means varying Xr​s and Xs​r simultaneously if r≠s) is:

∂ln​det(𝚿k)∂(𝚿k)r​s =(2−δr​s)​(𝚿k−1)r​s=(2−δr​s)​(𝚺k)r​s
∂(𝒚n−𝝁k)⊤​𝚿k​(𝒚n−𝝁k)∂(𝚿k)r​s =(2−δr​s)​(𝒚n−𝝁k)r​(𝒚n−𝝁k)s

Therefore,

∂ln⁡Φk​(𝒚n;𝜶)∂(𝚿k)r​s=12​(2−δr​s)​(𝚺k)r​s−12​(2−δr​s)​(𝒚n−𝝁k)r​(𝒚n−𝝁k)s

So,

∂ℓ​(𝒚n;𝜶)∂(𝚿k)r​s=tk​(𝒚n;𝜶)​(2−δr​s)2​((𝚺k)r​s−(𝒚n−𝝁k)r​(𝒚n−𝝁k)s)

The score component at 𝜶∗ is:

S(𝚿k)r​s​(𝒚n;𝜶∗)=−tk​(𝒚n;𝜶∗)​(2−δr​s)2​((𝚺k∗)r​s−(𝒚n−𝝁k∗)r​(𝒚n−𝝁k∗)s).

∎

For the two lemmas below, we will use the following assumptions:

Assumption 7 (Identifiability & Smoothness).

The GMM density fclust​(𝐲;𝛂) is identifiable. ℓ1​(𝐲;𝛂). is three times continuously differentiable w.r.t. 𝛂 in an open ball 𝔹​(𝛂∗,r0) around 𝛂∗.

Assumption 8 (Compact Parameter Space & True Parameter Properties).

𝜶∗ is an interior point of a compact set 𝚯𝕍′⊂𝔹​(𝛂∗,r0). This implies:

  • •

    Mixing proportions: πk∗≥πmin>0 for all k=1,…,K, for some constant πmin∈(0,1/K].

  • •

    Means: ‖𝝁k∗‖2≤η<∞ for all k.

  • •

    Covariance and Precision Matrices: The covariance matrices 𝚺k∗ have eigenvalues λ​(𝚺k∗) such that 0<sm≤λ​(𝚺k∗)≤sM<∞. Consequently, for the precision matrices 𝚿k∗=(𝚺k∗)−1, their eigenvalues λ​(𝚿k∗) satisfy 0<1/sM≤λ​(𝚿k∗)≤1/sm<∞. We define σmin2=sm, σmax2=sM. And for precision matrices, θmin=1/sM, θmax=1/sm.

All 𝛂∈𝚯𝕍′ satisfy these bounds. Let L3 be an upper bound on the norm of the third derivative tensor of ℓ1​(𝐲;𝛂) for 𝛂∈𝚯𝕍′, such that 𝔼​[L3​(𝐲)]<∞.

Assumption 9 (Data Distribution).

Observations 𝐲1,…,𝐲N are i.i.d. from fclust​(𝐲;𝛂∗). For each 𝐲n, there exists a latent class variable zn∈{1,…,K} with P​(zn=k)=πk∗, such that 𝐲n|zn=k∼𝒩​(𝛍k∗,𝚺k∗).

Assumption 10 (Bounded Posteriors).

For 𝛂∈𝚯𝕍′, the posterior probabilities tk​(𝐲;𝛂)=πk​ϕ​(𝐲|𝛍k,𝚺k)∑j=1Kπj​ϕ​(𝐲|𝛍j,𝚺j) satisfy 0<tk​(𝐲;𝛂)≤1.

The empirical Hessian is 𝑯N​(𝜶)=∇2ℓN​(𝜶). The Fisher Information Matrix is 𝑰​(𝜶∗)=𝔼𝜶∗​[∇2ℓ​(𝒚;𝜶∗)]. Moreover, we have the following definition of the true support of penalized parameters:

Definition 3 (True Support of Penalized Parameters 𝕊0).

Let 𝛂∗=(𝛑∗,{𝛍k∗}k=1K,{𝚿k∗}k=1K) be the true GMM parameter vector. Consider the penalty

P​(𝜶)=λ​∑k=1K‖𝝁k‖1+ρ​∑k=1K‖𝚿k‖1.

We define the true support set 𝕊0 as the collection of indices corresponding to nonzero parameters in 𝛂∗ that are subject to penalization:

𝕊0=𝕊𝝁∗∪𝕊𝚿∗,s0=|𝕊0|=s𝝁∗+s𝚿∗.

Here:

  1. 1.

    Support of means

    𝕊𝝁∗={(k,d):1≤k≤K, 1≤d≤p,(𝝁k∗)d≠0},s𝝁∗=|𝕊𝝁∗|.
  2. 2.

    Support of precision matrices

    𝕊𝚿∗={(k,r,s):1≤k≤K, 1≤r,s≤p,(𝚿k∗)r​s≠0},s𝚿∗=|𝕊𝚿∗|.
Assumption 11 (Restricted Eigenvalue Condition).

For any parameter increment vector 𝚫 indexed compatibly with 𝛂, we 𝚫𝕊0 for its restriction to indices in 𝕊0 and 𝚫(𝕊0)c for its complement. For c0≥1, define the cone

𝒞​(c0,𝕊0)={𝚫:∥𝚫(𝕊0)c∥1≤c0​∥𝚫𝕊0∥1}.

Assume the FIM 𝐈​(𝛂∗) is positive definite, then there exists κI>0 such that, for all 𝚫∈𝒞​(c0,𝕊0),

𝚫⊤​𝑰​(𝜶∗)​𝚫≥κI​‖𝚫‖22.
Remark 1.

The constant c0 is chosen to match the cone where the estimation error 𝛂^−𝛂∗ lies. A sufficient choice is c0≥A0+1A0−1 when the regularization level satisfies λ≥A0​‖∇ℓN​(𝛂∗)‖∞ with A0>1.

Assumption 12 (Uniform Hessian Concentration).

et dα=(K−1)+K​D+K​D​(D+1)/2 denote the number of free parameters in 𝛂. There exist constants CH>0, c1H,c2H>0 and a radius δR>0 such that, for N≳s0​ln⁡dα, with probability at least 1−c1H​dα−c2H,

sup𝜶~∈𝔹​(𝜶∗,δR)∩𝚯𝕍′𝚫∈𝒞​(c0,𝕊0),‖𝚫‖2=1|𝚫⊤​(𝑯N​(𝜶~)−𝑰​(𝜶~))​𝚫|≤CH​s0​ln⁡dαN.

The constant CH depends on the bounds in Assumption 8 and on moment bounds for derivatives of ℓ. The radius can be taken as δR≍s0​ln⁡dα/N.

Lemma 2 (Gradient Bound).

Let 𝛂∗=(𝛑∗,𝛍1∗,…,𝛍K∗,𝚿1∗,…,𝚿K∗) be the true parameter vector of a K-component GMM with D-dimensional components. Suppose Assumptions 7-10 hold. Let dα be defined in Assumption 12. Then there exist absolute constants C1,C2>0 and Cg>0 such that, for any N≥1,

ℙ​(‖∇ℓN​(𝜶∗)‖∞≤Cg​ln⁡dαN)≥ 1−C1​dα−C2,

where ℓN​(𝛂)=1N​∑n=1Nℓ​(𝐲n;𝛂) is the empirical negative log-likelihood (recall ℓ​(𝐲;𝛂)=−ln⁡fclust​(𝐲;𝛂)). An explicit choice is Cg=2​(C2+1)/c0​νmax, where c0 is an absolute constant in Bernstein’s inequality (e.g., c0=1/2). This bound holds provided N≥CB​ln⁡dα, with

CB=2​(C2+1)​νmax2c0​(minj:bj≠0⁡(𝔼​[Sj​n2]/bj))2,

νmax2=maxj⁡𝔼​[Sj​n2], and bj is the sub-exponential scale of Sj​n (bj=0 for sub-Gaussian components).

Proof of Lemma 2.

By Assumptions 7 and 9 (regularity and correct specification),

𝔼𝜶∗​[Sj​(𝒚n;𝜶∗)]=0 for all ​j=1,…,dα

Let νj2=𝔼​[Sj​(𝒚n;𝜶∗)2] be the variance of the score evaluated at αj∗. We bound variances and tail parameters for each parameter block.

Mixing proportions Sπk​(yn;α∗). By Lemma 1,

Sπk​(𝒚n;𝜶∗)=tK​(𝒚n;𝜶∗)πK∗−tk​(𝒚n;𝜶∗)πk∗.

By Assumption 10, 0<tm​(𝒚n;𝜶∗)≤1, and by Assumption 8, πm∗≥πmin>0; hence |Sπk​(𝒚n;𝜶∗)|≤2/πmin. Thus Sπk is bounded sub-Gaussian with νπk2=𝔼​[Sπk2]≤(2/πmin)2 and bπk=0.

Mean parameters S(μk)d​(yn;α∗). From Lemma 1,

S(𝝁k)d​(𝒚n;𝜶∗)=−tk​(𝒚n;𝜶∗)​(𝚿k∗​(𝒚n−𝝁k∗))d.

Let Un​k​d:=(𝚿k∗​(𝒚n−𝝁k∗))d. Conditionally on zn=ℓ, we have 𝒚n=𝝁ℓ∗+εℓ with εℓ∼𝒩​(𝟎,𝚺ℓ∗). Then

Var​(Un​k​d∣zn=ℓ)=ed⊤​𝚿k∗​𝚺ℓ∗​𝚿k∗​ed≤‖𝚿k∗‖22​‖𝚺ℓ∗‖2≤θmax2​σmax2,

and

|𝔼[Un​k​d∣zn=ℓ]|=|ed⊤𝚿k∗(𝝁ℓ∗−𝝁k∗)|≤∥𝚿k∗∥2∥𝝁ℓ∗−𝝁k∗∥2≤θmax(2η).

Hence 𝔼​[Un​k​d2]≤θmax2​(σmax2+4​η2). Since |tk|≤1,

ν(𝝁k)d2=𝔼​[S(𝝁k)d2]≤𝔼​[Un​k​d2]≤θmax2​(σmax2+4​η2).

Moreover, Un​k​d is (conditionally) Gaussian and hence sub-Gaussian; with the bounded multiplier tk, S(𝝁k)d is sub-Gaussian, so b(𝝁k)d=0.

Precision entries S(𝚿k)r​s​(yn;α∗). From Lemma 1,

S(𝚿k)r​s​(𝒚n;𝜶∗)=−tk​(𝒚n;𝜶∗)​(2−δr​s)2​((𝚺k∗)r​s−(𝒚n−𝝁k∗)r​(𝒚n−𝝁k∗)s).

Let Vn​k​r​s:=(𝚺k∗)r​s−(𝒚n−𝝁k∗)r​(𝒚n−𝝁k∗)s. Each (𝒚n−𝝁k∗)r is sub-Gaussian with parameter ≲σmax+η, so the product (𝒚n−𝝁k∗)r​(𝒚n−𝝁k∗)s is sub-exponential; hence Vn​k​r​s is sub-exponential. Since |tk|≤1, S(𝚿k)r​s is sub-exponential. Thus there exist constants CΨ,ν,CΨ,b>0 such that

ν(𝚿k)r​s2=𝔼​[S(𝚿k)r​s2]≤CΨ,ν​(σmax2+η2)2,b(𝚿k)r​s≤CΨ,b​(σmax2+η2).

Union bound. Let νmax2=maxj⁡𝔼​[Sj​n2] and bmax=maxj⁡{bj:bj≠0}; by the above bounds these depend only on (πmin,η,sm,sM). By Bernstein’s inequality, for each coordinate j and any t>0,

ℙ​(|1N​∑n=1NSj​n|≥t)≤ 2​exp⁡(−c0​N​min⁡(t2νj2,tbj)),

with the convention that if bj=0 then min⁡(⋅)=t2/νj2. Set tN=Cg​ln⁡dαN and assume N≥CB​ln⁡dα so that tN≤minj:bj≠0⁡νj2/bj. Then for all j,

ℙ​(|1N​∑n=1NSj​n|≥tN)≤ 2​exp⁡(−c0​N​tN22​νj2)≤ 2​exp⁡(−c0​N​tN22​νmax2)=2​dα−c0​Cg22​νmax2.

A union bound over j=1,…,dα yields

ℙ​(‖∇ℓN​(𝜶∗)‖∞≥tN)≤ 2​dα 1−c0​Cg22​νmax2.

Choosing C1=2 and Cg so that c0​Cg22​νmax2=C2+1 gives the claim, i.e., Cg=2​(C2+1)/c0​νmax. Finally, the condition N≥CB​ln⁡dα is guaranteed by taking

CB=Cg2(minj:bj≠0⁡νj2/bj)2=2​(C2+1)​νmax2c0​(minj:bj≠0⁡(𝔼​[Sj​n2]/bj))2.

∎

Lemma 3 (Parameter Consistency for Penalized GMM Estimator).

Let 𝛂^ be any local minimizer of the penalized negative average log-likelihood

Q​(𝜶)=ℓN​(𝜶)+P​(𝜶),

where

ℓN​(𝜶)=−1N​∑n=1Nln⁡fclust​(𝒚n;𝜶),P​(𝜶)=λ​∑k=1K‖𝝁k‖1+ρ​∑k=1K‖𝚿k‖1.

Let 𝛂∗ be the true GMM parameter vector, and 𝕊0 be the true support of the penalized parameters in 𝛂∗, with sparsity s0=|𝕊0|. Let dα be the total number of free parameters in 𝛂 as defined in Assumptions 12.

Suppose Assumptions 7-12 hold and the gradient bound in Lemma 2 is satisfied. Choose the regularization parameters λ=ρ=λchosen, where

λchosen=A0​Cg​ln⁡dαN

for a constant A0>1 (e.g. A0=3). Then, provided N≳s0​ln⁡dα, with probability at least 1−C1grad​dα−C2grad−c1H​dα−c2H:

  1. 1.

    (L2-norm consistency):

    ‖𝜶^−𝜶∗‖2≤2​(A0+1)​CgκI​s0​ln⁡dαN.
  2. 2.

    (L1-norm consistency):

    ‖𝜶^−𝜶∗‖1≤4​A0​(A0+1)​Cg(A0−1)​κI​s0​ln⁡dαN.

Here Cg is the gradient bound constant from Lemma 2, and κI is the restricted eigenvalue constant from Assumption 11. The factor A0 is a user-chosen constant controlling the regularization strength.

Proof of Lemma 3.

Let 𝜶^ be a minimizer of the penalized negative log-likelihood Q​(𝜶)=ℓN​(𝜶)+P​(𝜶), where P​(𝜶)=λ​∑k‖𝝁k‖1+ρ​∑k‖𝚿k‖1 and ℓN​(𝜶)=1N​∑n=1Nℓ​(𝒚n;𝜶) with ℓ​(𝒚;𝜶)=−ln⁡fclust​(𝒚;𝜶). Since 𝜶^ minimizes Q, we have

ℓN​(𝜶^)−ℓN​(𝜶∗)≤P​(𝜶∗)−P​(𝜶^). (20)

Let 𝚫=𝜶^−𝜶∗ and assume 𝜶^∈𝚯𝕍′ so that 𝜶∗+𝚫∈𝚯𝕍′. A second-order Taylor expansion gives

ℓN​(𝜶^)−ℓN​(𝜶∗)=⟨∇ℓN​(𝜶∗),𝚫⟩+12​𝚫⊤​𝑯N​(𝜶∗+t0​𝚫)​𝚫

for some t0∈(0,1). Write 𝜶~=𝜶∗+t0​𝚫. We lower bound 12​𝚫⊤​𝑯N​(𝜶~)​𝚫 by decomposing

𝚫⊤​𝑯N​(𝜶~)​𝚫 =𝚫⊤​𝑰​(𝜶∗)​𝚫+𝚫⊤​(𝑰​(𝜶~)−𝑰​(𝜶∗))​𝚫+𝚫⊤​(𝑯N​(𝜶~)−𝑰​(𝜶~))​𝚫.

From this point, all bounds are stated for arbitrary 𝚫 in the cone 𝒞​(c0,𝕊0); the fact that the actual error 𝜶^−𝜶∗ lies in 𝒞​(c0,𝕊0) will be established later by Lemma 4.

By Assumption 11, for 𝚫∈𝒞​(c0,𝕊0),

𝚫⊤​𝑰​(𝜶∗)​𝚫≥κI​‖𝚫‖22.

By Assumption 7 and compactness (Assumption 8), 𝜶↦𝑰​(𝜶)=𝔼​[∇2ℓ​(𝒚;𝜶)] is Lipschitz on 𝔹​(𝜶∗,r0)∩𝚯𝕍′ in spectral norm: there exists LI>0 such that

‖𝑰​(𝜶~)−𝑰​(𝜶∗)‖2≤LI​‖𝜶~−𝜶∗‖2=LI​t0​‖𝚫‖2≤LI​‖𝚫‖2,

and therefore

|𝚫⊤​(𝑰​(𝜶~)−𝑰​(𝜶∗))​𝚫|≤‖𝑰​(𝜶~)−𝑰​(𝜶∗)‖2​‖𝚫‖22≤LI​‖𝚫‖23≤LI​δR​‖𝚫‖22,

provided ‖𝚫‖2≤δR. By Assumption 12, if ‖𝚫‖2≤δR and 𝚫∈𝒞​(c0,𝕊0), then with probability at least 1−c1H​dα−c2H,

|𝚫⊤​(𝑯N​(𝜶~)−𝑰​(𝜶~))​𝚫|≤CH​s0​ln⁡dαN​‖𝚫‖22.

Combining the three representations, for 𝚫∈𝒞​(c0,𝕊0) with ‖𝚫‖2≤δR,

12​𝚫⊤​𝑯N​(𝜶~)​𝚫≥12​(κI−LI​‖𝚫‖2−CH​s0​ln⁡dαN)​‖𝚫‖22.

Choose δR and N so that LI​δR≤κI/4 and CH​s0​ln⁡dαN≤κI/4 (e.g., N≳(CH2/κI2)​s0​ln⁡dα). Then, with the same probability,

12​𝚫⊤​𝑯N​(𝜶~)​𝚫≥κI4​‖𝚫‖22.

Hence we obtain the restricted strong convexity (RSC) inequality on the cone:

ℓN​(𝜶∗+𝚫)−ℓN​(𝜶∗)−⟨∇ℓN​(𝜶∗),𝚫⟩≥κL2​‖𝚫‖22,for all ​𝚫∈𝒞​(c0,𝕊0),‖𝚫‖2≤δR, (21)

where κL=κI/2>0.

Stochastic term. By Lemma 2, if N≥CB​ln⁡dα, then with probability at least 1−C1grad​dα−C2grad,

‖∇ℓN​(𝜶∗)‖∞≤Cg​ln⁡dαN.

By Hölder’s inequality (⟨x,y⟩<=∥x∥∞​∥y∥1),

|⟨∇ℓN​(𝜶∗),𝚫⟩|≤‖∇ℓN​(𝜶∗)‖∞​‖𝚫‖1≤Cg​ln⁡dαN​‖𝚫‖1. (22)

Penalty difference. Let 𝚫𝝁k=𝝁^k−𝝁k∗ and 𝚫𝚿k=𝚿^k−𝚿k∗, and denote by 𝕊𝝁∗ and 𝕊𝚿∗ the true supports (Definition 3). Using the L1-triangle inequality on supports: for any vector 𝒂,𝒃 and support 𝕊 of 𝒂: ‖𝒂‖1−‖𝒃‖1≤‖𝒂𝕊−𝒃𝕊‖1−‖𝒂(𝕊)c−𝒃(𝕊)c‖1+2​‖𝒂(𝕊)c‖1. Since 𝒂(𝕊)c=𝟎:

‖𝝁k∗‖1−‖𝝁^k‖1≤‖𝚫𝝁k,𝕊𝝁∗‖1−‖𝚫𝝁k,(𝕊𝝁∗)c‖1,‖𝚿k∗‖1−‖𝚿^k‖1≤‖𝚫𝚿k,𝕊𝚿∗‖1−‖𝚫𝚿k,(𝕊𝚿∗)c‖1.

Summing over k and writing 𝕊0=𝕊𝝁∗∪𝕊𝚿∗,

P​(𝜶∗)−P​(𝜶^) ≤λ​(‖𝚫𝝁,𝕊0‖1−‖𝚫𝝁,(𝕊0)c‖1)+ρ​(‖𝚫𝚿,𝕊0‖1−‖𝚫𝚿,(𝕊0)c‖1). (23)

Combining. With probability at least 1−C1grad​dα−C2grad−c1H​dα−c2H, combining Equation 21, Equation 22, and Equation 23 with Equation 20 yields, for all 𝚫∈𝒞​(c0,𝕊0) with ‖𝚫‖2≤δR,

κL2​‖𝚫‖22≤Cg​ln⁡dαN​‖𝚫‖1+λ​(‖𝚫𝝁,𝕊0‖1−‖𝚫𝝁,(𝕊0)c‖1)+ρ​(‖𝚫𝚿,𝕊0‖1−‖𝚫𝚿,(𝕊0)c‖1). (24)

Invoking the inequality Equation 29 from Lemma 4 (the cone lemma, proved via KKT and the gradient bound, thus independent of RSC),

κL2​‖𝚫‖22≤(1+1A0)​∑j∈𝕊0λj​|Δj|−(1−1A0)​∑j∈(𝕊0)cλj​|Δj|.

Since the second term is non-positive (as A0>1 and λj​|Δj|≥0), we can drop it to get an upper bound:

κL2​‖𝚫‖22≤(1+1A0)​∑j∈𝕊0λj​|Δj|.

Taking the common regularization level λchosen=λ=ρ=A0​Cg​ln⁡dαN with A0>1 and note that ∑j∈𝕊0λj​|Δj|=λchosen​‖𝚫𝕊0‖1, we obtain

κL2​‖𝚫‖22≤(1+1A0)​λchosen​‖𝚫𝕊0‖1.

Using the Cauchy-Schwarz inequality ‖𝚫𝕊0‖1≤s0​‖𝚫𝕊0‖2, and since ‖𝚫𝕊0‖2≤‖𝚫‖2:

κL2​‖𝚫‖22≤λchosen​(1+1A0)​s0​‖𝚫‖2.

If ‖𝚫‖2≠0, we can divide by ‖𝚫‖2:

‖𝚫‖2≤2​λchosenκL​(1+1A0)​s0.

Substituting λchosen=A0​Cg​ln⁡dαN:

‖𝜶^−𝜶∗‖2 ≤2​A0​CgκL​(1+1A0)​s0​ln⁡dαN
=2​(A0+1)​CgκL​s0​ln⁡dαN.

This proves the L2-rate with constant CL​2=2​(A0+1)​CgκL.

For the L1-rate, Lemma 4 yields the cone bound ‖𝚫(𝕊0)c‖1≤Ccone​‖𝚫𝕊0‖1 with Ccone=A0+1A0−1 (decomposability via [44]). Hence, assuming all λj for penalized components are λchosen):

‖𝚫‖1 =‖𝚫𝕊0‖1+‖𝚫(𝕊0)c‖1
≤‖𝚫𝕊0‖1+Ccone​‖𝚫𝕊0‖1=(1+Ccone)​‖𝚫𝕊0‖1
≤(1+Ccone)​s0​‖𝚫𝕊0‖2(by Cauchy-Schwarz)
≤(1+Ccone)​s0​‖𝚫‖2.

Substituting Ccone=A0+1A0−1:

1+Ccone=1+A0+1A0−1=A0−1+A0+1A0−1=2​A0A0−1.

So,

‖𝚫‖1≤(2​A0A0−1)​s0​‖𝚫‖2.

Now substitute the bound for ‖𝚫‖2:

‖𝜶^−𝜶∗‖1 ≤(2​A0A0−1)​s0​(2​(A0+1)​CgκL​s0​ln⁡dαN)
=4​A0​(A0+1)​Cg(A0−1)​κL​s0​ln⁡dαN.
=CL​1​(A0,Cg,κL)​s0​ln⁡dαN

where CL​1​(A0,Cg,κL)=4​A0​(A0+1)​Cg(A0−1)​κL. This concludes the proof of parameter consistency in L1 and L2 norms. ∎

Remark 2.

The radius δR governing the RSC inequality Equation 21 is critical and is typically of the same order as the target statistical error. If ‖𝚫‖2 exceeds δR, the Lipschitz remainder LI​‖𝚫‖23 can dominate, and/or the uniform Hessian concentration in Assumption 12 may fail on such a large neighborhood of 𝛂∗. In general M-estimation analyses (e.g., [44]), one often proves an RSC with a tolerance term,

ℓN​(𝜶∗+𝚫)−ℓN​(𝜶∗)−⟨∇ℓN​(𝜶∗),𝚫⟩≥κL2​‖𝚫‖22−τL​ln⁡dαN​‖𝚫‖12,

valid on a larger set. Under Assumption 12, we obtain sufficiently strong control to dispense with this tolerance and derive the cleaner quadratic curvature bound Equation 21. If Assumption 12 were weakened (e.g., only yielding a bound of the form CH​ln⁡dαN​‖𝚫‖1s0​‖𝚫‖2 on 𝚫⊤​(𝐇N−𝐈)​𝚫 over the cone), then a tolerance term proportional to ‖𝚫‖12 would naturally appear in the RSC.

Lemma 4 (Cone Condition for L1-Penalized M-Estimators).

Let 𝛂^ be any local minimizer of Q​(𝛂)=ℓN​(𝛂)+P​(𝛂), where ℓN​(𝛂) is a differentiable loss function and the penalty P​(𝛂) is a sum of component-wise L1 penalties:

P​(𝜶)=∑j=1dαλj​|αj|.

Let 𝕊0 be the true support of 𝛂∗. Let 𝚫=𝛂^−𝛂∗. Suppose the regularization parameters are chosen such that for some constant A0>1:

λj≥A0​|[∇ℓN​(𝜶∗)]j|for all ​j∈(𝕊0)c. (25)

and λj of similar order for j∈𝕊0. For simplicity, we often set λj=λchosen for all penalized components, where λchosen≥A0​‖∇ℓN​(𝛂∗)‖∞. If the RSC condition from Equation 21 holds for 𝚫:

ℓN​(𝜶∗+𝚫)−ℓN​(𝜶∗)−⟨∇ℓN​(𝜶∗),𝚫⟩≥κL2​‖𝚫‖22,

then, for A0>1, the error vector 𝚫 satisfies the cone condition:

∑j∈(𝕊0)cλj​|Δj|≤A0+1A0−1​∑j∈𝕊0λj​|Δj|. (26)

If all λj for penalized components are equal to λchosen, this simplifies to:

‖𝚫(𝕊0)c‖pen,1≤A0+1A0−1​‖𝚫𝕊0‖pen,1,

where ∥⋅∥pen,1 refers to the L1 norm over the components that are actually penalized. If all components were penalized with the same λ0, this would be ‖𝚫(𝕊0)c‖1≤A0+1A0−1​‖𝚫𝕊0‖1.

Remark 3.

In our case, λj=λ for mean components μk​d, and λj=ρ for off-diagonal precision components (Ψk)r​s, and λj=0 for parameters not penalized like proportions or diagonal precision elements if they are not penalized towards a specific value. Moreover, the Cone Condition ensures that the error vector is primarily concentrated on the true support. The key is that the regularization parameter for the "noise" variables (off-support) must be sufficiently larger than the corresponding component of the score vector, allowing the penalty to effectively shrink noise components.

Proof of Lemma 4.

Let 𝒫⊂{1,…,dα} denote the index set of penalized coordinates and write ‖𝒙‖pen,1=∑j∈𝒫|xj|. In what follows, sums and L1-norms are taken over 𝒫. Assume the tuning condition holds on all penalized coordinates:

λj≥A0​|[∇ℓN​(𝜶∗)]j|∀j∈𝒫,A0>1, (27)

equivalently λchosen≥A0​‖∇ℓN​(𝜶∗)‖∞,pen when a common level is used. 222If some coordinates are unpenalized, they are excluded from 𝒫 and do not enter the cone inequality.

Since 𝜶^ is a local minimizer of Q​(𝜶), it satisfies the basic optimality inequality

ℓN​(𝜶^)−ℓN​(𝜶∗)≤P​(𝜶∗)−P​(𝜶^).

Invoking the restricted curvature bound (RSC) Equation 21 for 𝚫=𝜶^−𝜶∗,

ℓN​(𝜶^)−ℓN​(𝜶∗)≥⟨∇ℓN​(𝜶∗),𝚫⟩+κL2​‖𝚫‖22.

Note that, only the nonnegativity of the quadratic remainder is needed for the cone; the explicit κL>0 is used later for L2-rates. Therefore,

⟨∇ℓN​(𝜶∗),𝚫⟩+κL2​‖𝚫‖22≤P​(𝜶∗)−P​(𝜶^). (28)

For the penalty difference, write P​(𝜶)=∑j∈𝒫λj​|αj| and let 𝕊0⊆𝒫 denote the true support of the penalized coordinates. Using the standard L1 support inequality with 𝒂=𝜶∗ and 𝒃=𝜶^=𝜶∗+𝚫,

‖𝒂‖1−‖𝒃‖1=‖𝒂𝕊0‖1−‖𝒃𝕊0‖1−‖𝒃(𝕊0)c‖1≤‖𝒂𝕊0−𝒃𝕊0‖1−‖𝒃(𝕊0)c‖1,

and since (𝜶∗)(𝕊0)c=𝟎 on 𝒫, we obtain

P​(𝜶∗)−P​(𝜶^)≤∑j∈𝕊0λj​|Δj|−∑j∈(𝕊0)cλj​|Δj|.

Substitute this into Equation 28 and bound the linear term by the triangle inequality over 𝒫:

⟨∇ℓN​(𝜶∗),𝚫⟩≤∑j∈𝒫|[∇ℓN​(𝜶∗)]j|​|Δj|≤∑j∈𝒫λjA0​|Δj|,

where the last step uses Equation 27. We obtain

κL2​‖𝚫‖22≤∑j∈𝕊0(λj+λjA0)​|Δj|−∑j∈(𝕊0)c(λj−λjA0)​|Δj|.

Equivalently,

κL2​‖𝚫‖22≤(1+1A0)​∑j∈𝕊0λj​|Δj|−(1−1A0)​∑j∈(𝕊0)cλj​|Δj|. (29)

Since the left-hand side is nonnegative, it follows that

(1−1A0)​∑j∈(𝕊0)cλj​|Δj|≤(1+1A0)​∑j∈𝕊0λj​|Δj|,

and for A0>1 this yields the cone inequality

∑j∈(𝕊0)cλj​|Δj|≤A0+1A0−1​∑j∈𝕊0λj​|Δj|.

If all penalized coordinates share a common level λchosen, this becomes ‖𝚫(𝕊0)c‖pen,1≤A0+1A0−1​‖𝚫𝕊0‖pen,1, and if every coordinate is penalized equally it reduces to ‖𝚫(𝕊0)c‖1≤A0+1A0−1​‖𝚫𝕊0‖1. ∎

With the cone condition in Lemma 4 established on the penalized coordinates, we now consolidate the standing assumptions for the ranking-consistency analysis (Lemma 5) and for Theorem 3. The goal is to avoid duplication, align the SRUW assumptions used earlier with the GMM penalized framework here, and make explicit exactly which conditions are invoked downstream.

Remark 4 (Alignment with SRUW assumptions).

Assumptions 4-6 were introduced for the SRUW model in the MNARz setting. In the present penalized GMM analysis for Theorem 3:

  • •

    Model uniqueness (SRUW-4). We condition on a fixed mixture structure (number of components K and dimension D) and require identifiability of the GMM density; this role is played by Assumption 7 below. We do not require the SRUW tuple (K0,m0,r0,l0,𝕍0) explicitly here because no model selection over SRUW structures is performed in Theorem 3.

  • •

    Compactness and interiority (SRUW-5-6). These correspond directly to Assumption 8 below (compact parameter subset 𝚯𝕍′ and 𝜶∗ interior), together with the eigenvalue and boundedness constraints therein. Thus SRUW-5-6 are subsumed by Assumption 8.

Hence, for the purposes of Lemma 5 and Theorem 3, it suffices to work with Assumptions 7-12 below; SRUW-4-6 need not be re-stated.

To proceed rigorously towards Lemma  5 and, ultimately, Theorem 3, we now restate the full set of standing assumptions-consolidating those aligned with the 𝕊​ℝ​𝕌​𝕎 framework and those introduced earlier for penalized likelihood analysis-into a unified assumption block.

Assumption 13 (Standing Assumptions for proving Theorem 3).

The following holds:

  1. 1.

    Identifiability & Smoothness (Assumption 7): the GMM density fclust​(𝒚;𝜶) is identifiable and ℓ​(𝒚;𝜶)=−ln⁡fclust​(𝒚;𝜶) is three times continuously differentiable in an open ball 𝔹​(𝜶∗,r0).

  2. 2.

    Compactness & Bounds (Assumption 8): 𝜶∗ is an interior point of a compact 𝚯𝕍′⊂𝔹​(𝜶∗,r0); mixture weights, means, and covariance/precision eigenvalues satisfy the stated uniform bounds (with πmin,η,σmin2,σmax2,θmin,θmax).

  3. 3.

    Data-generating mechanism (Assumption 9): 𝒚1,…,𝒚N i.i.d. from fclust​(⋅;𝜶∗) with latent zn∼Mult​(𝝅∗) and 𝒚n|zn=k∼𝒩​(𝝁k∗,𝚺k∗).

  4. 4.

    Posterior responsibilities (Assumption 10): tk​(𝒚;𝜶)∈(0,1] for all 𝜶∈𝚯𝕍′.

  5. 5.

    Fisher RE on a cone (Assumption 11): restricted eigenvalue condition for 𝑰​(𝜶∗) on 𝒞​(c0,𝕊0) with constant κI>0.

  6. 6.

    Uniform Hessian concentration (Assumption 12): for radius δR≍s0​ln⁡dα/N and N≳s0​ln⁡dα,

    sup𝜶~∈𝔹​(𝜶∗,δR)∩𝚯𝕍′sup𝚫∈𝒞​(c0,𝕊0),‖𝚫‖2=1|𝚫⊤​(𝑯N​(𝜶~)−𝑰​(𝜶~))​𝚫|≤CH​s0​ln⁡dαN

    with probability at least 1−c1H​dα−c2H.

  7. 7.

    Penalty on penalized coordinates. On the penalized index set 𝒫, choose λj so that

    λj≥A0​‖∇ℓN​(𝜶∗)‖∞,penfor some ​A0>1,

    e.g. a common level λchosen=A0​Cg​ln⁡dαN with Cg from Lemma 2. Unpenalized coordinates are excluded from 𝒫 and from all ∥⋅∥pen,1 norms.

Assumption 14 (Working high-probability event for ranking analysis).

Let ℰgrad={‖∇ℓN​(𝛂∗)‖∞,pen≤Cg​ln⁡dα/N} be the event in Lemma 2, and ℰRSC the event on which Equation 21 holds with curvature κL>0 over {𝚫∈𝒞​(c0,𝕊0):‖𝚫‖2≤δR}. Under Assumption 13 and for N large enough, both events hold with probability at least 1−δN, where δN=C1grad​dα−C2grad+c1H​dα−c2H. In the proof of Lemma 5, we condition on ℰgrad∩ℰRSC.

Definition 4 (Quantities for Ranking Consistency).

Let 𝛂∗ be the true GMM parameters, and 𝛂^​(λ) be the estimator for a given regularization level λ (with ρ either tied to λ or fixed).

  • •

    Noise level.

    λnoise=Cg​ln⁡dαN,

    the uniform bound on the score at 𝜶∗ from Lemma 2.

  • •

    Effective curvature for means.

    Hk​jeff​(𝜶∗)=𝔼𝜶∗​[tk​(𝒚;𝜶∗)​(𝚿k∗)j​j]=(𝚿k∗)j​j​𝔼𝜶∗​[tk​(𝒚;𝜶∗)],

    where the expectation is with respect to

    𝒚∼fclust​(⋅;𝜶∗).

    By Assumption 8, 𝔼𝜶∗​[tk​(𝒚;𝜶∗)]=πk∗≥πmin and (𝚿k∗)j​j≥θmin (since 𝚿k∗≻0 and ej⊤​𝚿k∗​ej≥λmin​(𝚿k∗)), hence

    Hk​jeff​(𝜶∗)≥πmin​θmin> 0.
  • •

    KKT remainder for mean coordinates.

    Rk​j​(𝜶^​(λ),𝜶∗,λ)=∂μk​jℓN​(𝜶^​(λ))−∂μk​jℓN​(𝜶∗)−Hk​jeff​(𝜶∗)​(μ^k​j​(λ)−μk​j∗).

    We assume that, with high probability and for λ in the regime λ≍(ln⁡dα)/N,

    |Rk​j​(𝜶^​(λ),𝜶∗,λ)|≤EKKT_rem​(λ),

    where EKKT_rem​(λ) depends on ‖𝜶^​(λ)−𝜶∗‖ (cf. Lemma 3) and on bounds for Hessian/third derivatives (from Assumptions 7-12). In particular, when ‖𝚫‖2 is of order s0​ln⁡dα/N, we take

    EKKT_rem​(λ)≤Crem​λnoise

    for some constant Crem≥0 over the relevant range of λ.

  • •

    Noise-screening threshold.

    ΛN∗=(1+Crem+ϵN)​λnoise,

    with a small margin ϵN>0.

  • •

    Signal-preservation threshold. For a coordinate j and some k0 attaining (or meeting) a signal condition for |μk0​j∗|,

    ΛS∗​(j;λ)=Hk0​jeff​(𝜶∗)​|μk0​j∗|−(1+ϵS)​λnoise−(1+ϵS)​EKKT_rem​(λ),

    where ϵS>0 is a fixed margin. Choosing λ so that ΛS∗​(j;λ)>0 ensures the signal at (k0,j) persists (i.e., μ^k0​j​(λ)≠0) under the KKT inequalities.

Assumption 15 (True Relevant/Irrelevant Variables).

Let 𝕊𝛍∗ denote the set of indices j such that at least one component mean has a nonzero j-th entry, i.e., j∈𝕊𝛍∗ iff ∃k∈{1,…,K} with μk​j∗≠0. For j∉𝕊𝛍∗, we have μk​j∗=0 for all k.

Assumption 16 (Minimum Signal Strength).

For each j∈𝕊𝛍∗, there exists k0∈{1,…,K} such that, for all λ up to some detection level λdetect​_​upper,

Hk0​jeff​(𝜶∗)​|μk0​j∗|≥λ+(1+ϵS)​λnoise+(1+ϵS)​EKKT​_​rem​(λ), (30)

with ϵS>0 fixed. In particular, Equation 30 implies ΛS∗​(j;λ)>λ. Moreover, we assume the separation

minj∈𝕊𝝁∗⁡λj†>ΛN∗,where ​λj†:=inf{λ′>0:λ′=ΛS∗​(j;λ′)}, (31)

i.e., the smallest regularization at which the j-th signal would vanish strictly exceeds the noise-screening threshold ΛN∗.

Assumption 17 (Regularization Grid).

The grid 𝒢λ={λ(1),…,λ(MG)} covers a neighborhood of the noise threshold ΛN∗ and extends up to

minj∈𝕊𝝁∗⁡λj†with ​λj†​ as in Equation 31.

Assume ρ is either tied to λ (e.g., ρ≍λ) or fixed so that its effect is absorbed by constants in the bounds.

Remark 5.

The remainder Rk​j​(𝛂^,𝛂∗,λ) in Definition 4 captures the deviation of the mean-coordinate KKT equation from its linearized form. It aggregates: (i) off-diagonal Fisher blocks acting on other coordinates of 𝚫, (ii) empirical-population curvature fluctuations (𝐇N−𝐈), and (iii) population curvature drift 𝐈​(𝛂~)−𝐈​(𝛂∗) from evaluating at 𝛂~. Under the standing assumptions (smoothness/compactness: Assumptions 7-8, Hessian concentration: Assumption 12) and the L2-error bound from Lemma 3, one obtains the bound

|Rk​j​(𝜶^​(λ),𝜶∗,λ)|≤EKKT​_​rem​(λ)≤Crem​λnoise

on the high-probability event ℰgrad∩ℰRSC, provided the sparsity/sample-size regime ensures CH​s0​ln⁡dαN+LI​‖𝛂^​(λ)−𝛂∗‖2≲1 (e.g., s0≲N/ln⁡dα). A crude but sufficient bound for individual coordinates is ‖𝚫‖∞≤‖𝚫‖2=O​(s0​ln⁡dα/N) by Lemma 3. Assumption 16 then ensures the effective signal Hk0​jeff​(𝛂∗)​|μk0​j∗| dominates both the regularization level and the stochastic/remainder terms, yielding persistence of true signals and suppression of noise along the regularization path.

Lemma 5 (Ranking Consistency for Mean Parameters).

Under Assumptions 13, 15, 16, and 17, and on the high-probability event ℰgrad∩ℰRSC from Assumption 14, let the grid 𝒢λ contain points distributed across [0,λgrid_max] with

λgrid_max>minj∈𝕊𝝁∗⁡λj†,λj†:=inf{λ′>0:λ′=ΛS∗​(j;λ′)}.

Then, with probability at least

Prank≥ 1−MG⋅D⋅δN,where ​δN:=C1grad​dα−C2grad+c1H​dα−c2H,

MG=|𝒢λ|, and D denotes the number of penalized mean coordinates per variable (e.g., D=K), the following hold:

  1. 1.

    Relevant Variables. For each j∈𝕊𝝁∗,

    𝒪K​(j)≥Nsignal​(j):=#​{λ∈𝒢λ:λ<ΛS∗​(j;λ)}.

    Let ηR:=minj∈𝕊𝝁∗⁡Nsignal​(j). By the separation Equation 31,

    ηR>#​{λ∈𝒢λ:λ≤ΛN∗}.
  2. 2.

    Irrelevant Variables. For each j∉𝕊𝝁∗,

    𝒪K​(j)≤Nnoise:=#​{λ∈𝒢λ:λ≤ΛN∗}.

Consequently, ηR>Nnoise, yielding a strict separation in ranking scores between relevant and irrelevant variables with probability at least Prank.

Proof of Lemma 5.

Work on the high-probability event

ℰrank:=ℰgrad∩ℰRSC

(on which the uniform gradient bound |Sk​j∗|≤λnoise and the KKT-remainder bound |Rk​j​(𝜶^​(λ),𝜶∗,λ)|≤EKKT_rem​(λ) hold simultaneously for all penalized mean coordinates (k,j) and all λ∈𝒢λ). By Assumptions 13-17 and the union bound,

ℙ​(ℰrank)≥ 1−MG⋅D⋅δN,δN:=C1grad​dα−C2grad+c1H​dα−c2H.

The ranking score for variable j is

𝒪K(j)=∑λ∈𝒢λℐ(∃k∈{1,…,K}:μ^k​j(λ)≠0).

KKT expansion. For a penalized mean coordinate (k,j), the KKT condition reads

0=∇μk​jℓN​(𝜶^​(λ))+λ​ξ^k​j,ξ^k​j∈{{sign⁡(μ^k​j​(λ))},μ^k​j​(λ)≠0,[−1,1],μ^k​j​(λ)=0.

A first-order expansion at 𝜶∗ yields

∇μk​jℓN​(𝜶^​(λ))=[∇μk​jℓN​(𝜶∗)]⏟=⁣:Sk​j∗+Hk​jeff​(𝜶∗)​(μ^k​j​(λ)−μk​j∗)+Rk​j​(𝜶^​(λ),𝜶∗,λ),

hence

λ​ξ^k​j=−Sk​j∗−Hk​jeff​(𝜶∗)​(μ^k​j​(λ)−μk​j∗)−Rk​j​(𝜶^​(λ),𝜶∗,λ). (32)

Relevant variables (j∈𝕊μ∗). Fix j∈𝕊𝝁∗. By Assumption 16 there exists k0 with μk0​j∗≠0 satisfying the signal condition. Assume, for contradiction, that at a given λ∈𝒢λ we have μ^k0​j​(λ)=0. Then ξ^k0​j∈[−1,1] and μ^k0​j​(λ)−μk0​j∗=−μk0​j∗. Applying Equation 32 and taking absolute values,

λ≥|Hk0​jeff​(𝜶∗)​μk0​j∗−Sk0​j∗−Rk0​j|≥Hk0​jeff​(𝜶∗)​|μk0​j∗|−|Sk0​j∗|−|Rk0​j|.

On ℰrank,

λ≥Hk0​jeff​(𝜶∗)​|μk0​j∗|−λnoise−EKKT_rem​(λ).

Therefore, whenever

λ​<Hk0​jeff​(𝜶∗)|​μk0​j∗|−λnoise−EKKT_rem​(λ),

we must have μ^k0​j​(λ)≠0. Since ΛS∗​(j;λ)=Hk0​jeff​(𝜶∗)​|μk0​j∗|−(1+ϵS)​λnoise−(1+ϵS)​EKKT_rem​(λ) is a stricter threshold,

λ<ΛS∗​(j;λ)⟹μ^k0​j​(λ)≠0.

Hence the count of grid points for which variable j is (at least in one component) active satisfies

𝒪K​(j)≥Nsignal​(j):=#​{λ∈𝒢λ:λ<ΛS∗​(j;λ)}.

Taking the minimum over j∈𝕊𝝁∗ gives ηR:=minj∈𝕊𝝁∗⁡Nsignal​(j).

Irrelevant variables (j∉𝕊μ∗). Here μk​j∗=0 for all k. For a given λ, the zero solution μ^k​j​(λ)=0 is KKT-feasible iff

|∇μk​jℓN​(𝜶^​(λ))|≤λ.

Using the expansion with μk​j∗=0 and μ^k​j​(λ)=0,

|Sk​j∗+Rk​j​(𝜶^​(λ),𝜶∗,λ)|≤λ.

On ℰrank, this is guaranteed whenever

λ≥λnoise+EKKT_rem​(λ).

By definition of the noise threshold ΛN∗=(1+ϵN)​λnoise+(1+ϵN)​supλ′≤ΛN∗EKKT_rem​(λ′) and the monotone/worst-case domination in its definition, any λ>ΛN∗ satisfies λ≥λnoise+EKKT_rem​(λ). Therefore, for each j∉𝕊𝝁∗ and all λ>ΛN∗, all coordinates μ^k​j​(λ) equal zero, so the indicator ℐ(∃k:μ^k​j(λ)≠0)=0. Consequently,

𝒪K​(j)≤Nnoise:=#​{λ∈𝒢λ:λ≤ΛN∗}.

Separation and probability. By the separation Assumption Equation 31,

ηR>Nnoise,

yielding a strict gap between the scores of relevant and irrelevant variables on ℰrank. Finally, by the union bound over at most MG grid points and D penalized mean coordinates per variable, we obtain

ℙ​(the above conclusions hold for all ​j)≥ 1−MG⋅D⋅δN,

which is Prank in the statement. This completes the proof. ∎

Theorem 4 (Selection Consistency of the Two-Step SRUW Procedure).

Assume all assumptions in Lemma 3 and Lemma 5, together with Theorem 2 for the final (K,m,r,ℓ) choice, hold. Let s0=|𝕊0| and MN​R=p−s0 be the number of non-𝕊0 variables. Moreover, suppose the following assumptions hold:

  1. (a)

    False negative for 𝕊0. For any true relevant variable j∈𝕊0 and any intermediate set 𝕊^cur⊂𝕊0 with j∉𝕊^cur, the probability of incorrectly rejecting j from 𝕊 by the BICdiff criterion is uniformly bounded:

    sup𝕊^cur⊂𝕊0ℙ​(BICdiff​(j∣𝕊^cur)≤0)≤pS​(N),

    where s0⋅pS​(N)≤ϵS,F​N​(N) and ϵS,F​N​(N)=oN​(1).

  2. (b)

    False positive for non-𝕊0. For any true non-relevant variable j∉𝕊0 (i.e., j∈𝕌0∪𝕎0), given that 𝕊0 has been correctly identified (i.e., 𝕊^cur=𝕊0), the probability of BICdiff​(j∣𝕊0)>0 is bounded:

    ℙ​(BICdiff​(j∣𝕊0)>0)≤pN​(N),

    where pN​(N)→0 as N→∞. For BIC, typically pN​(N)=𝒪​(N−γB) for some γB>0 (e.g., γB≥Δ​νmin/2 where Δ​νmin≥1 is the minimum parameter-penalty gap).

  3. (c)

    Consistency of BIC-penalized regressions for ℝ^​[j∣𝕊0].:

    • •

      For j∈𝕎0, ℙ​(ℝ^​[j∣𝕊0]=∅)≥1−preg​(N).

    • •

      For j∈𝕌0, ℙ​(ℝ^​[j∣𝕊0]=ℝ0​(j)≠∅)≥1−preg​(N).

    where (w0+u0)​preg​(N)=oN​(1).

Let the stopping-rule parameter for selecting 𝕊^ (and 𝕎^) be c, understood so that once 𝕊^cur=𝕊0 holds, each j∉𝕊0 is tested at most c times by the BICdiff(⋅∣𝕊0) screening before termination.

For a desired tolerance ϵS,F​P>0 on the probability of including any false positive in 𝕊^, it suffices to choose c so that

MN​R​(1−(1−pN​(N))c)≤ϵS,F​P⟺c≤ln⁡(1−ϵS,F​P/MN​R)ln⁡(1−pN​(N)). (33)

For small pN​(N), this is well-approximated by

c≲ϵS,F​PMN​R​pN​(N).

In particular, if pN​(N) decays polynomially in N (e.g., N−γB with γB≥1/2) and MN​R​pN​(N)→0, then any fixed cfixed≥1 (e.g., cfixed=1 or 3) ensures

ℙ​(any false positive in ​𝕊^)≤MN​R​(1−(1−pN​(N))cfixed)≤MN​R​cfixed​pN​(N)→ 0.

Then, under (a)-(c) and with c chosen according to Equation 33 (or with c=cfixed such that MN​R​cfixed​pN​(N)→0), the two-step SRUW procedure recovers the true model structure (K0,m0,r0,ℓ0,𝕍0) with probability ℙ​(Success)→1 as N→∞.

Proof.

We show each stage succeeds with high probability and then combine them.

Conditioning on ranking accuracy. By Lemma 5, with probability at least 1−Pfail,rank_total, the ranking event ℰrank holds: all s0 relevant variables in 𝕊0 appear before any pure-noise 𝕎0 variables, and 𝕌0 appear after 𝕊0. We condition on ℰrank for the remainder of the argument.

No false negatives in 𝕊^. Fix j∈𝕊0 and any 𝕊^cur⊂𝕊0 not containing j. By Assumption (a),

ℙ​(BICdiff​(j∣𝕊^cur)≤0)≤pS​(N).

By a union bound over the s0 true variables,

ℙ(𝕊0⊈𝕊^|ℰrank)≤s0pS(N)=:ϵS,F​N(N)→0.

Let ℰS-noFN be the event 𝕊0⊆𝕊^. Then ℙ​(ℰS-noFN∣ℰrank)≥1−ϵS,F​N​(N).

No false positives in 𝕊^. On ℰrank∩ℰS-noFN we have 𝕊^cur=𝕊0. For any j∉𝕊0, Assumption (b) gives

ℙ(BICdiff(j∣𝕊0)>0)≤pN(N)=:q.

Let ℰS-noFP be the event that no non-relevant variable is ever added to 𝕊^ after 𝕊0 is reached.

(A) Distribution-free bound. For any fixed set 𝒥 of candidates examined of any size,

ℙ(∃j∈𝒥:BICdiff(j∣𝕊0)>0)≤|𝒥|q

by a union bound. In particular,

ℙ​(ℰS-noFPc|ℰrank∩ℰS-noFN)≤MN​R​q,

which tends to 0 if MN​R​pN​(N)→0. This bound requires no independence and is always valid.

(B) Sharper bound under c-run termination. Assume the screening step proceeds through non-relevant candidates until it observes c consecutive rejections, and that the c decisions in such a run are (asymptotically) independent or, more generally, satisfy

ℙ​(⋂t=1c{BICdiff​(jt∣𝕊0)≤0})≥(1−q)c

for any c distinct non-relevant candidates (j1,…,jc). Then the probability that a given c-block fails (i.e., contains at least one FP) is

ℙ​(block error)≤ 1−(1−q)c.

Partition the MN​R non-relevant indices into ⌊MN​R/c⌋ disjoint blocks of size c (discarding a remainder if needed). By a union bound over blocks,

ℙ​(ℰS-noFPc|ℰrank∩ℰS-noFN)≤MN​Rc​(1−(1−q)c)≤MN​R​(1−(1−q)c).

Therefore, enforcing

MN​R​(1−(1−pN​(N))c)≤ϵS,F​P⟺c≤ln⁡(1−ϵS,F​P/MN​R)ln⁡(1−pN​(N))

yields ℙ​(ℰS-noFPc∣ℰrank∩ℰS-noFN)≤ϵS,F​P. For small pN​(N), 1−(1−pN​(N))c∼c​pN​(N), so c≲ϵS,F​P/(MN​R​pN​(N)).

Combining (A)-(B): whenever the independence/decoupling condition for runs holds, we may use the sharper (1−(1−q)c) design in the theorem statement; otherwise, the distribution-free guarantee MN​R​pN​(N) is valid (and implies the sharper one whenever c is fixed and MN​R​pN​(N)→0).

Put together 𝕊^=𝕊0. Let ℰS∗:=ℰS-noFN∩ℰS-noFP. Then, with either (A) or (B),

ℙ​(𝕊^=𝕊0|ℰrank)≥ 1−ϵS,F​N​(N)−ϵS,F​P​(N)→ 1,

for ϵS,F​P​(N)=MN​R​pN​(N) in (A), or ϵS,F​P​(N)=MN​R​(1−(1−pN​(N))c) in (B).

Consistency of 𝕎^,𝕌^,ℝ^. Given 𝕊^=𝕊0, the reverse scan uses BIC-penalized regressions to decide 𝕎^ and ℝ^. By Assumption (c),

ℙ​(𝕎^=𝕎0|ℰrank∩ℰS∗)≥ 1−(w0+u0)​preg​(N)−ϵW,F​P​(N)→ 1,

and

ℙ​(ℝ^=ℝ0|ℰrank∩ℰS∗∩{𝕎^=𝕎0})≥ 1−u0​preg​(N)→ 1.

Final SRUW choice via BIC. By Theorem 2, the final selection of (K,m,r,ℓ) is consistent, i.e., ℙ​(Pfinal_choice)→1.

Multiplying the stage-wise success probabilities,

ℙ​(Overall Success)≥ℙ​(ℰrank)⋅ℙ​(ℰS∗∣ℰrank) ⋅ℙ​(𝕎^=𝕎0∣⋯)
⋅ℙ(ℝ^=ℝ0∣⋯)⋅ℙ(Pfinal_choice)→ 1.

The slowest decaying term among Pfail,rank_total,ϵS,F​N​(N),ϵS,F​P​(N),(w0+u0)​preg​(N),u0​preg​(N) governs the overall rate. In particular, under the sharper design (B) with c chosen as above,

ϵS,F​P​(N)=MN​R​(1−(1−pN​(N))c)≤ϵS,F​P,

while the distribution-free design (A) yields ϵS,F​P​(N)=MN​R​pN​(N)→0 whenever MN​R​pN​(N)→0. ∎

Theorem 5 (Equivalence of MNARz-SRUW and MAR on Augmented Data for SRUW).

Consider an observation (𝐲n,𝐜n) where 𝐲n=(𝐲nS,𝐲nU,𝐲nW) and 𝐜n is its missingness pattern. Let the complete-data likelihood for 𝐲n given cluster zn​k=1 under the SRUW model (K,m,r,ℓ,𝕍) be:

fSRUW​(𝒚n∣zn​k=1;θk)=fk​(𝒚nS;𝜶k)​freg​(𝒚nU∣𝒚nR;𝜽reg)​findep​(𝒚nW;𝜽indep).

Assume an MNARz-SRUW mechanism for missing data, where the probability of the missingness pattern 𝐜n depends only on the cluster membership zn​k:

f​(𝒄n∣𝒚n,zn​k=1;ψk) =f​(𝒄n∣zn​k=1;ψk) (34)
=∏d∈S∪U∪Wρk​dcn​d​(1−ρk​d)1−cn​d,

where ρk​d∈(0,1) are components of the missingness parameter ψk (thus f​(𝐜n∣zn​k=1;ψk) does not depend on 𝐲n). The observed data under this MNARz-SRUW model is (𝐲no,𝐜n), and its likelihood is:

LMNARz-SRUW​(𝒚no,𝒄n;𝜽,𝝍)=∫∑k=1Kπk​fSRUW​(𝒚n∣zn​k=1;θk)​f​(𝒄n∣zn​k=1;ψk)​d​𝒚nm. (35)

Now consider the augmented observed vector 𝐲~no=(𝐲no,𝐜n). Assume 𝐲~no is i.i.d. from a mixture model in which the conditional law of the missing part 𝐲nm is MAR with respect to 𝐲~no (i.e., the missingness mechanism does not depend on 𝐲nm given (𝐲no,𝐜n)). Define the **augmented observed-data** likelihood as

f~MAR(𝒚~no;𝜽,𝝍)=∑k=1Kπk(∫fSRUW(𝒚no,𝒚nm∣zn​k=1;θk)d𝒚nm)f(𝒄n∣zn​k=1;ψk), (36)

where fSRUW(⋅∣zn​k=1;θk) is the same component density as above.

Then, for fixed parameters (𝛉,𝛙), the observed-data likelihood under MNARz-SRUW for (𝐲no,𝐜n) is identical to the likelihood of the augmented observation 𝐲~no=(𝐲no,𝐜n) under the MAR interpretation:

LMNARz-SRUW​(𝒚no,𝒄n;𝜽,𝝍)=f~MAR​(𝒚~no;𝜽,𝝍).
Proof of Theorem 5.

Starting from Equation 35, by the MNARz assumption f​(𝒄n∣𝒚n,zn​k=1;ψk)=f​(𝒄n∣zn​k=1;ψk) is independent of 𝒚nm, hence

LMNARz-SRUW​(𝒚no,𝒄n;𝜽,𝝍) =∫∑k=1KπkfSRUW(𝒚no,𝒚nm∣zn​k=1;θk)f(𝒄n∣zn​k=1;ψk)d𝒚nm
=∑k=1Kπkf(𝒄n∣zn​k=1;ψk)(∫fSRUW(𝒚no,𝒚nm∣zn​k=1;θk)d𝒚nm)
=∑k=1Kπk​fk​(𝒚no∣zn​k=1;θk)​f​(𝒄n∣zn​k=1;ψk),

where we set fk(𝒚no∣zn​k=1;θk):=∫fSRUW(𝒚no,𝒚nm∣zn​k=1;θk)d𝒚nm. Comparing with Equation 36 yields

LMNARz-SRUW​(𝒚no,𝒄n;𝜽,𝝍)=f~MAR​(𝒚~no;𝜽,𝝍),

as claimed. ∎

By Theorem 5, the observed-data likelihood under MNARz-SRUW equals the augmented observed-data likelihood under MAR for every (𝜽,𝝍) and each observation; hence, the full-sample likelihoods coincide pointwise in (𝜽,𝝍). The following corollary shows the resulting estimator-level equivalence.

Corollary 1 (Estimator equivalence under augmentation).

Under the conditions of Theorem 5, for any fixed (𝛉,𝛙) and for each observation,

LMNARz-SRUW​(𝒚no,𝒄n;𝜽,𝝍)=f~MAR​(𝒚~no;𝜽,𝝍).

Hence the full-sample observed-data likelihoods coincide, and therefore the sets of maximum likelihood estimators for (𝛉,𝛙) under MNARz-SRUW based on (𝐲o,𝐜) and under MAR on the augmented data 𝐲~o=(𝐲o,𝐜) are identical (up to label switching). Moreover, if the same priors on (𝛉,𝛙) are used in both formulations, the resulting Bayesian posteriors coincide.

Proof of Corollary 1.

By Theorem 5, for every observation n and every fixed (𝜽,𝝍),

LMNARz-SRUW​(𝒚no,𝒄n;𝜽,𝝍)=f~MAR​(𝒚~no;𝜽,𝝍),𝒚~no=(𝒚no,𝒄n).

Let 𝒟={(𝒚no,𝒄n)}n=1N denote the sample, assumed i.i.d. under the mixture model as in the theorem. The full-sample observed-data likelihoods are then

LMNARz-SRUW​(𝒟;𝜽,𝝍) =∏n=1NLMNARz-SRUW​(𝒚no,𝒄n;𝜽,𝝍)
L~MAR​(𝒟;𝜽,𝝍) =∏n=1Nf~MAR​(𝒚~no;𝜽,𝝍).

By pointwise equality for each factor, we obtain the full-sample equality

LMNARz-SRUW​(𝒟;𝜽,𝝍)=L~MAR​(𝒟;𝜽,𝝍)for all ​(𝜽,𝝍).

Consequently, the sets of maximum likelihood estimators coincide:

arg⁡max(𝜽,𝝍)⁡LMNARz-SRUW​(𝒟;𝜽,𝝍)=arg⁡max(𝜽,𝝍)⁡L~MAR​(𝒟;𝜽,𝝍),

up to the usual permutations of mixture component labels (label switching), since both likelihoods are invariant under relabelings.

For the Bayesian statement, let Π be a common prior on (𝜽,𝝍) admitting a density π​(𝜽,𝝍) with respect to a common dominating measure. The posteriors are

πMNARz​(𝜽,𝝍∣𝒟) ∝π​(𝜽,𝝍)​LMNARz-SRUW​(𝒟;𝜽,𝝍)
πMAR​(𝜽,𝝍∣𝒟) ∝π​(𝜽,𝝍)​L~MAR​(𝒟;𝜽,𝝍).

Since the likelihoods coincide pointwise in (𝜽,𝝍), the unnormalized posteriors agree, hence their normalizing constants (integrals over the same parameter space) are also equal. Therefore

πMNARz​(𝜽,𝝍∣𝒟)=πMAR​(𝜽,𝝍∣𝒟),

again, modulo label switching. ∎

Remark 6.

In the MAR interpretation of the augmented data 𝐲~no, the components 𝐲no have density fk​(𝐲no∣zn​k=1;θk) in cluster k, and the components 𝐜n (which are fully observed within 𝐲~no) have density f​(𝐜n∣zn​k=1;ψk) in cluster k. The MAR assumption applies to 𝐲nm, meaning its missingness mechanism, given 𝐲no and 𝐜n, does not depend on 𝐲nm itself. The likelihood of 𝐲~no is formed by integrating out 𝐲nm from f​(𝐲n,𝐜n∣zn​k=1), which leads to the product fk​(𝐲no∣zn​k=1;θk)​f​(𝐜n∣zn​k=1;ψk) for each cluster k.

Appendix E Detailed EM Algorithms for the Two-Stage Procedure

This section gives complete EM derivations for both stages. In Stage A (ranking), we run a penalized GMM on the imputed and standardized data 𝒀¯=std​(𝒀~), where 𝒀~ is the single-imputed matrix defined in Algorithm 1 which is used for ranking only. In Stage B (role assignment), we fit the unpenalized SRUW model under MNARz by exploiting the equivalence to MAR on the augmented data (𝒀,𝑪). Throughout, 𝚿k=𝚺k−1, the SRUW partition is V=(𝕊,ℝ,𝕌,𝕎), and tn​k denotes EM responsibilities. We write ℓ​(⋅) for (penalized) log-likelihoods and Q​(⋅;⋅) for EM Q-functions.

Stage A (ranking): Penalized EM for adaptive-regularized GMM

Penalized objective.

Given 𝒚¯1,…,𝒚¯N∈ℝD, we maximize the penalized observed-data log-likelihood

ℓpen​(𝜶)=∑n=1Nlog⁡[∑k=1Kπk​ϕ​(𝒚¯n∣𝝁k,𝚺k)]−λ​∑k=1K‖𝝁k‖1−ρ​∑k=1K∑i≠j𝑷k,i​j​|𝚿k,i​j|. (37)

Assumption: For each EM run at a fixed (λ,ρ), the weights 𝑷k are computed from the warm start and then held fixed.333If 𝑷k is updated during EM, one must add a majorization step to preserve monotonicity.

At iteration t, with responsibilities tn​k(t), the penalized Q-function is

Qpen​(𝜶;𝜶(t−1)) =∑n=1N∑k=1Ktn​k(t)​{log⁡πk−12​log⁡|𝚺k|−12​(𝒚¯n−𝝁k)⊤​𝚿k​(𝒚¯n−𝝁k)}
−λ​∑k=1K‖𝝁k‖1−ρ​∑k=1K∑i≠j𝑷k,i​j​|𝚿k,i​j|. (38)
E-step.
tn​k(t)=πk(t−1)​ϕ​(𝒚¯n∣𝝁k(t−1),𝚺k(t−1))∑ℓ=1Kπℓ(t−1)​ϕ​(𝒚¯n∣𝝁ℓ(t−1),𝚺ℓ(t−1)),nk(t)=∑n=1Ntn​k(t). (39)
M-step: mixing weights.

πk(t)=nk(t)/N.

M-step: means 𝝁k with ℓ1 penalty.

Given 𝚿k(t−1), the subproblem in 𝝁k is convex. Let

𝒎¯k(t):=1nk(t)​∑n=1Ntn​k(t)​𝒚¯n.

Then the gradient of the smooth part is

∇𝝁k[12​∑n=1Ntn​k(t)​(𝒚¯n−𝝁k)⊤​𝚿k(t−1)​(𝒚¯n−𝝁k)]=nk(t)​𝚿k(t−1)​(𝝁k−𝒎¯k(t)).

The KKT optimality for coordinate j is

nk(t)​(𝚿k(t−1))j​j​μk​j+nk(t)​∑v≠j(𝚿k(t−1))j​v​μk​v−nk(t)​(𝚿k(t−1)​𝒎¯k(t))j∈λ​∂|μk​j|.

A coordinate-descent update is

μk​j←1nk(t)​(𝚿k(t−1))j​j​𝒮λ​(nk(t)​(𝚿k(t−1)​𝒎¯k(t))j−nk(t)​∑v≠j(𝚿k(t−1))j​v​μk​v),

where 𝒮λ​(u)=sign​(u)​max⁡{|u|−λ,0}.

M-step: precisions 𝚿k via weighted graphical lasso.

With 𝝁k(t) fixed, define the responsibility-weighted covariance

𝑺k(t)=1nk(t)​∑n=1Ntn​k(t)​(𝒚¯n−𝝁k(t))​(𝒚¯n−𝝁k(t))⊤.

Then

𝚿k(t)∈arg⁡min𝚿≻0⁡{−log​det𝚿+tr​(𝑺k(t)​𝚿)+2​ρnk(t)​∑i≠j𝑷k,i​j​|𝚿i​j|},

with diagonals unpenalized; standard glasso solvers apply.

Stage B (role assignment): EM for SRUW under MNARz

Observed likelihood via augmentation.

Let 𝒟MAR​∪˙​𝒟MNAR=[D]. Under MNARz,

ℓ​(Θ;𝒀,𝑪)=∑n=1Nlog⁡[∑k=1Kπk​fk,MARo​(𝒚no;𝜶k,𝝃)​fcMNARz​(𝒄n,MNAR;𝝍k)],

is a standard mixture on the augmented observation (𝒚no,𝒄n,MNAR).

Complete-data log-likelihood and Q-function.

With Θ=(𝝅,{𝜶k},𝝃,{𝝍k}) and latent {zn​k},

Q​(Θ;Θ(t−1))=∑n=1N∑k=1Ktn​k(t)​{log⁡πk+𝔼​[log⁡fk​(𝒀n;𝜶k,𝝃)∣𝒚no,zn​k=1]+log⁡fcMNARz​(𝒄n,MNAR;𝝍k)}.
E-step.
tn​k(t)=πk(t−1)​fk,MARo​(𝒚no;𝜶k(t−1),𝝃(t−1))​fcMNARz​(𝒄n,MNAR;𝝍k(t−1))∑ℓ=1Kπℓ(t−1)​fℓ,MARo​(𝒚no;𝜶ℓ(t−1),𝝃(t−1))​fcMNARz​(𝒄n,MNAR;𝝍ℓ(t−1)).
M-step: πk and MNARz parameters.
πk(t)=1N​∑n=1Ntn​k(t),ρ^k(t)=∑n=1Ntn​k(t)​∑d∈𝒟MNARcn​d∑n=1Ntn​k(t)​|𝒟MNAR|,𝝍k(t)=log⁡ρ^k(t)1−ρ^k(t).
M-step: SRUW data model updates.

Write

fk​(𝒚;𝜶k,𝝃)=fclust​(𝒚𝕊;𝝁k𝕊,𝚺k𝕊)​freg​(𝒚𝕌∣𝒚ℝ;𝒂,𝜷,𝛀)​findep​(𝒚𝕎;𝜸,𝚪).

Let 𝔼k[⋅∣𝒚no] denote Gaussian conditional expectations under component k. Then the sufficient statistics are the mixture-weighted moments

𝔼[⋅∣𝒚no]=∑k=1Ktn​k(t)𝔼k[⋅∣𝒚no,zn​k=1].

Cluster block 𝕊: for each k, compute

𝒚¯k𝕊,(t)=1nk(t)​∑n=1Ntn​k(t)​𝔼k​[𝒚n𝕊∣𝒚no],𝑺k𝕊,(t)=1nk(t)​∑n=1Ntn​k(t)​𝔼k​[𝒚n𝕊​(𝒚n𝕊)⊤∣𝒚no],

and set 𝝁k𝕊,(t)=𝒚¯k𝕊,(t) and 𝚺k𝕊,(t)=𝑺k𝕊,(t)−𝒚¯k𝕊,(t)​(𝒚¯k𝕊,(t))⊤.

Regression block 𝕌∣ℝ: let 𝑿n=[ 1(𝒚nℝ)⊤]. Form the mixture-weighted global moments

𝑨=∑n=1N∑k=1Ktn​k(t)​𝔼k​[𝑿n⊤​𝑿n∣𝒚no],𝑩=∑n=1N∑k=1Ktn​k(t)​𝔼k​[𝑿n⊤​𝒚n𝕌∣𝒚no],

then [𝒂(t)𝜷(t)]=𝑨−1​𝑩, and

𝛀(t)=1N​∑n=1N∑k=1Ktn​k(t)​𝔼k​[(𝒚n𝕌−𝒂(t)−𝒚nℝ​𝜷(t))​(𝒚n𝕌−𝒂(t)−𝒚nℝ​𝜷(t))⊤|𝒚no].

Independent block 𝕎:

𝜸(t)=1N​∑n=1N∑k=1Ktn​k(t)​𝔼k​[𝒚n𝕎∣𝒚no],𝚪(t)=1N​∑n=1N∑k=1Ktn​k(t)​𝔼k​[𝒚n𝕎​(𝒚n𝕎)⊤∣𝒚no]−𝜸(t)​(𝜸(t))⊤.
Numerical stability.

In Stage A, we use warm starts and pathwise (λ,ρ); diagonals unpenalized and enforce 𝚿k≻0. In Stage B, if some nk(t) is tiny, we add a small ridge to 𝚺k𝕊 or merge/discard components per BIC.

Appendix F Additional Experiments

F.1 Metric Details

We outline some metrics that we use to evaluate the clustering accuracy and imputation error, which are summarized below:

Table 3: Summary of evaluation metrics for clustering and imputation quality
Name Formula Range
Imputation error
NRMSE 1N​∑i=1N(yi−y^i)2σ [0,1]
WNRMSE ∑c=1C(NRMSEc×wc)∑c=1Cwc [0,1]
Clustering similarity
ARI RI−𝔼​[RI]max⁡(RI)−𝔼​[RI] [0,1]
Composite
CIIE α×(1−NRMSE)+𝜷×Similarity Score [0,1]

F.2 Assessing Integrated MNARz Handling

We conduct an additional experiment to showcase the performance of handling MNAR pattern within the variable selection framework of which the setups are similar to [62]. Datasets were simulated with n=100 observations, K=3 true clusters with proportions 𝝅=(0.5,0.25,0.25)), and D=6,9 variables. True cluster memberships 𝐙 were drawn according to 𝝅. The complete data 𝒀 was generated as Yn​d=∑k′=1KZn​k′​δk′​d+ϵn​d, where ϵn​d∼𝒩​(0,1), and 𝜹 defined cluster-specific mean shifts with signal strength τ=2.31 with δ11=δ14=δ22=δ25=δ33=δ36=τ, others zero). For a specific MNAR scenario, the class-specific intercept component ψkz was 0, and variable-specific slopes ψjy were (1.45,0.2,−3,1.45,0.2,−3), resulting in P​(Mi​j=1|Yi​j,Zi​k=1)=Φ​(ψjy​Yi​j).

In Figure 4 and 5, we plot the boxplots of ARI and NRMSE over 20 replications. We observe that MNAR-based approaches attain better results over methods not designed for MNAR patterns. These methods deliver competitive performance even when the true missingness mechanism is more complex (MNARy, MNARyz); thereby, supporting the conclusion in [62]. Furthermore, our framework demonstrates even slightly higher ARI and lower NRMSE in some cases compared to standard MNARz.

Refer to caption
Figure 4: Boxplot of the ARI obtained over 20 replications of simulated data. The theoretical ARIs are represented by a red dashed line.
Refer to caption
Figure 5: Boxplot of the NRMSE obtained over 20 replications of simulated data

F.3 Sensitivity on the choice of c

As discussed in [13] and further investigated through Theorem 4, the choice of the hyperparameter c plays a crucial role in the stepwise construction of the relevant variable set 𝕊^, balancing the risk of stopping the selection process too early with the risk of incorrectly including irrelevant variables. Our theoretical work, particularly Theorem 4 and the formulation for selecting c in Equation 33, suggests that c should ideally be determined by considering the probability of a single incorrect inclusion (pN​(N)), the number of non-relevant variables (MN​R), and a desired tolerance for overall false positives in 𝕊^ (ϵS,F​P). The experimental results presented here provide practical insights into this interplay and generally support our theoretical conclusions.

As illustrated in Figures 6 and 7, the two datasets (Dataset 1: D=7,s0=3; Dataset 2: D=14,s0=2) show different optimal ranges for c. This aligns with the theoretical expectation that c is not a universal constant but rather interacts with dataset-specific factors, which are encapsulated by pN​(N) and MN​R. For Dataset 1, the excellent performance with c=2 (achieving 100% correct selection of relevant variables) points to a very low pN​(N). This suggests that after the three true relevant variables are found, the remaining four non-𝕊0 variables consistently produce BICdiff≤0, making a small stopping value like c=2 both effective and safe. The continued strong performance at c=7 (90% correct) further indicates that pN​(N) is small enough that even scanning all MN​R=4 non-relevant variables rarely leads to a false inclusion. This situation is consistent with a theoretical scenario where pN​(N)≪1/MN​R, minimizing the risk of false positives regardless of c within a reasonable range.

In contrast, Dataset 2, which has more non-𝕊0 variables (MN​R=12), demonstrates a clearer trade-off. Good performance is observed for c=2 (80%) and c=3 (90%), implying that pN​(N) is still relatively small. However, the notable decline in performance when c=7 (only 5% correct 𝕊0 selection) empirically validates a key theoretical concern: if c is set too high and pN​(N) is not sufficiently close to zero, the algorithm examines more non-𝕊0 variables while awaiting c consecutive negative BICdiff values. Each additional variable examined increases the cumulative chance of a Type I error (a non-𝕊0 variable having BICdiff>0 by chance). If the likelihood of obtaining c correct rejections in a row, (1−pN​(N))c, diminishes significantly as more variables are processed, false inclusions become more probable. The poor result for c=7 in Dataset 2 suggests its pN​(N) value makes it unlikely to achieve seven consecutive correct rejections before a false positive occurs among the MN​R=12 candidates. This aligns with the theoretical relationship c≈ln⁡(MN​R/ϵS,F​P)/pN​(N): for a given ϵS,F​P, a higher pN​(N) or a larger MN​R would generally favor a smaller c to maintain that error tolerance, or a large fixed c might lead to a poorer effective ϵS,F​P. The finding that c=3 is effective for Dataset 2 is consistent with the heuristic used in prior work [13], suggesting it offers a practical compromise when pN​(N) is small but non-negligible, and MN​R is moderate. Furthermore, the relative stability of cluster number selection across different values of c suggests that c’s main influence is on the selection of variables for 𝕊, as predicted by theory, with more indirect effects on the overall model choice, primarily if 𝕊^ is significantly misidentified.

Refer to caption
Figure 6: Proportions of choosing correct number clusters and relevant variables for simulated dataset in Section 5 under varying c. We run the experiment over 50 replications
Refer to caption
Figure 7: Clustering performance of simulated dataset in Section 5 under varying c. We run the experiment over 50 replications

F.4 Computational Times and Scalability

Complexity Analysis

We count arithmetic operations up to absolute constants (big-Oh). Throughout this part we denote: N∈ℕ (samples), D∈ℕ (variables), K∈ℕ (mixture components). Let MEM be a uniform upper bound on the per-fit number of EM iterations until convergence (Assumption A2 below). In Stage A (ranking), the M-step contains a graphical-lasso solve per component with at most Mglasso outer iterations. Covariance inverses and log-determinants are computed by Cholesky factorizations.

Assumptions. We work under the following explicit conditions.

  • •

    A1. The ground-truth SRUW partition has Deff:=|𝕊true|+|𝕌true|≪D.

  • •

    A2. Each EM fit terminates in at most MEM iterations. For population/regularized EM, this is justified by established convergence rates: either sublinear convergence MEM=𝒪p​(1/N) [31] or geometric convergence MEM=𝒪​(log⁡(1/ϵ)) to achieve ϵ-accuracy [76]. Our analysis uses the more conservative bounded iteration assumption for clarity.

  • •

    A3. With probability 1−o​(1), the Stage-A ranking lists all Deff informative variables before any purely irrelevant variables (proved in Theorem 3 of the paper).

  • •

    A4. Let Cglasso​(d) denote the arithmetic cost of one graphical-lasso solve of size d×d. In general, Cglasso​(d)=Θ​(Mglasso​d3). Under connected-component decomposition with largest block size smax (as in [73, 41]), Cglasso​(d)=Θ​(Mglasso​∑cpc3)≤Θ​(Mglasso​d​smax2).

Per-iteration costs for a d-variate GMM. One EM iteration for a K-component Gaussian mixture in dimension d has:

  • •

    E-step: Responsibilities tn​k require evaluating Gaussian log-densities for all (n,k). With per-component 𝚺k−1 and log​det𝚺k fixed within the iteration, each evaluation uses a matrix-vector product and a quadratic form, costing Θ​(d2). Hence Θ​(N​K​d2).

  • •

    M-step (means, weights): Weighted sums are Θ​(N​K​d).

  • •

    M-step (covariances): Weighted second moments yield Θ​(N​K​d2). The per-component matrix factorization/inversion is Θ​(d3), hence Θ​(K​d3).

Therefore one EM iteration in dimension d costs

Θ​(N​K​d2)+Θ​(K​d3),

and one EM fit (up to convergence) costs

𝒞EM​(d)=Θ​(MEM​(N​K​d2+K​d3)). (40)

The classical SRUW selection starts from all D variables and eliminates one at a time. At step j (j=D,D−1,…,2), to remove one variable it evaluates j candidates; each evaluation requires an EM fit in dimension j−1. Using Equation 40, the total cost is

∑j=2Dj⋅𝒞EM​(j−1)=Θ​(MEM​∑j=2Dj​(N​K​(j−1)2+K​(j−1)3)).

Using the polynomial sums:

∑j=1Dj3=D2​(D+1)24=Θ​(D4)

and

∑j=1Dj4=D​(D+1)​(2​D+1)​(3​D2+3​D−1)30=Θ​(D5)

, we obtain the tight bound

𝒞stepwise=Θ​(MEM​(N​K​D4+K​D5)). (41)

The D5 term renders this approach impractical beyond moderate D.

Stage A (Ranking). For each (λ,ρ) on a grid of size Mgrid, we run a penalized EM. Per iteration: the E-step remains Θ​(N​K​D2); the M-step adds K graphical-lasso solves. By Assumption A4,

per iter cost=Θ​(N​K​D2)+Θ​(K​Cglasso​(D)).

Therefore the total Stage-A cost is

𝒞rank=Θ​(Mgrid​MEM​(N​K​D2+K​Cglasso​(D))). (42)

Two explicit regimes follow immediately from A4:

(Dense/general) Cglasso​(D)=Θ​(Mglasso​D3)
⇒𝒞rank=Θ​(Mgrid​MEM​(N​K​D2+K​Mglasso​D3)). (43)
(Connected components of size ​smax​) Cglasso​(D)=Θ​(Mglasso​D​smax2)
⇒𝒞rank=Θ​(Mgrid​MEM​(N​K​D2+K​Mglasso​D​smax2)). (44)

Stage B (Role assignment). A single forward/backward pass evaluates a constant number of SRUW-MNARz fits per newly considered variable. Under A3 the pass stops after Deff additions; hence the Stage-B complexity is

𝒞role=Θ​(MEM​(N​K​Deff2+K​Deff3)). (45)

Since Deff≪D (A1), 𝒞role is strictly lower order than 𝒞rank and Stage A dominates. Combining Equation 43-Equation 44 and Equation 45, the two-stage complexity is

𝒞two-stage=𝒞rank+𝒞role=Θ​(Mgrid​MEM​(N​K​D2+K​Cglasso​(D)))+o​(𝒞rank). (46)

From Equation 41 and Equation 46, the speedup factor satisfies

𝒞stepwise𝒞two-stage=Ω​(N​K​D4+K​D5Mgrid​(N​K​D2+K​Cglasso​(D))).

Two explicit lower bounds follow.

  • •

    If N≥Mglasso​D (E-step dominates Stage A): using Equation 43,

    𝒞stepwise𝒞two-stage=Ω​(D2Mgrid).
  • •

    If N<Mglasso​D (glasso dominates Stage A): using Equation 43,

    𝒞stepwise𝒞two-stage=Ω​(D2Mgrid​Mglasso).

Under connected-component sparsity with block size smax Equation 44 the second case strengthens to

𝒞stepwise𝒞two-stage=Ω​(D2Mgrid⋅min⁡{1,NMglasso​smax2}).

In all regimes the speedup is Ω​(D2/Mgrid), and strictly larger when sparsity (small smax) is present.

Empirical Validation

We report wall-clock times (seconds) for representative scenarios, confirming the predicted polynomial speedups:

Scenario N D K SelvarMNARz (s) Clustvarsel (s) Speedup
Varying D 750 15 4 12.2 190 ∼15×
Varying D 750 21 4 14.4 640 ∼44×
Varying D 750 27 4 15.5 2054 ∼132×
Varying N 1000 20 4 19.0 838 ∼44×
Varying K 750 20 12 9.77 1242 ∼127×

These measurements are consistent with the theory: the two-stage method scales like D2 in the dense case (and better under sparsity), whereas backward stepwise scales as D4-D5.

Dependency of Computational Times under Degree of Missingness

We analyze the computational dependency of our framework on the proportion of missing data, revealing an advantageous property: runtime decreases with increasing missing rates due to the efficient handling of incomplete data patterns in the EM algorithm.

The key insight stems from the E-step computational complexity in Stage B, which operates directly on the observed data patterns. Let Xn denote the number of observed entries in observation 𝒚n. The E-step cost for Gaussian mixture models scales as:

𝒞E-step=Θ​(∑n=1N∑k=1KXn2)=Θ​(K​N​𝔼​[X2])

For different missingness mechanisms:

  • •

    MCAR at rate r: X∼Binom​(D,1−r), yielding

    𝔼​[X2]=Var​(X)+𝔼​[X]2=D​(1−r)​r+D2​(1−r)2

    This produces a quadratic decrease in E-step work as r↑1.

  • •

    MAR/MNARz: The same complexity bound applies, with Xn representing the random count of observed coordinates per observation. Each per-record Gaussian density evaluation scales quadratically with observed entries, maintaining the Θ​(K​N​𝔼​[X2]) complexity.

Stage A employs efficient single imputation and remains largely insensitive to missing rates, while Stage B inherits the beneficial 𝔼​[X2] scaling. Consequently, total runtime decreases monotonically with increasing missingness rates.

We empirically verified this property by fixing (N,D,K)=(1000,14,4) and systematically increasing the missing data rate under two mixed missingness scenarios. The results, summarized in Table 4 and Figure 8, confirm the theoretical predictions.

Table 4: Runtime (seconds) with increasing missing rate under mixed mechanisms.
Missing Mechanism 20% 30% 50% 80%
Mixed (MAR+MNAR) 54.1 46.2 32.9 23.9
Mixed (MCAR+MAR+MNAR) 34.9 23.8 25.5 22.8
Refer to caption
Figure 8: Runtime as the proportion of missing data increases.

Three key observations emerge:

  1. 1.

    Runtime consistently decreases with higher missing rates, as predicted by the 𝔼​[X2] scaling. The MNARz block processes smaller observed patterns during E-step calculations, reducing computational burden for conditional expectations and complete-data sufficient statistics.

  2. 2.

    The slight variation at 30%-50% missingness in the three-mechanism scenario likely stems from random allocation of missing positions. When missing data concentrates in clustering variables (𝕊), the number of EM iterations for convergence may vary, but the overall inverse relationship persists.

  3. 3.

    This property provides significant practical benefits. Users can expect faster processing on datasets with higher missing rates-a valuable characteristic for real-world applications where extensive missing data is common.

The combination of mathematical analysis and empirical results demonstrates that our framework not only handles high missing rates effectively but also becomes more computationally efficient as missing data increases, making it particularly suitable for challenging real-world datasets with substantial missingness.

F.5 More on Simulated Dataset

Performance under Model Misspecification. Previous simulations focused on “pure” MAR and MNAR scenarios to clarify the benefits of our model when combined with the MNARz mechanism. However, real-world missingness mechanisms are often mixed. To assess robustness, we conducted additional simulations with a dataset of D=14 variables, where the true clustering variable set is 𝕊⋆={1,2} (remaining variables play roles of ℝ,𝕌,𝕎 according to scenario 8 in section 5. We generated missing data using a mixture of mechanisms:

  • •

    MCAR: Some variables have missing values completely at random

  • •

    MAR: Missingness depends only on observed components yno

  • •

    MNARy: Missingness depends directly on unobserved values 𝒚m (actively violating MNARz assumption)

The MNARy mechanism (see [62] for precise formulation) specifically tests our model’s robustness to misspecification, as it depends on the unobserved values rather than cluster assignments. We compare our method against several baselines, including a novel two-step approach designed to isolate the benefits of joint modeling:

  1. 1.

    MNARz + SelvarMix: Estimate GMM-MNARz model and impute missing data, then run SelvarMix [13] on the imputed dataset

  2. 2.

    Multiple Imputation variants: Using both missRanger and gcimputeR (cite here) with Mclust

  3. 3.

    VarSelLCM: A competing joint modeling approach

Table 5: Experiment with mixed MAR+MNARy (True set {1,2})
Method ARI NRMSE Relevant Variables Selected
SelvarMNARz (ours) 0.808 0.176 1, 2 (correct)
VarSelLCM 0.779 – 1-11 (extra variables)
missRanger + Mclust 0.799 0.190 –
gcimputeR + Mclust 0.774 0.413 –
MNARz + SelvarMix 0.328 0.070 1 (missing variable)
Table 6: Experiment with mixed MCAR+MAR+MNARy (True set {1,2})
Method ARI NRMSE Relevant Variables Selected
SelvarMNARz (ours) 0.808 0.181 1, 2 (correct)
VarSelLCM 0.772 – 1-11 (extra variables)
missRanger + Mclust 0.783 0.198 –
gcimputeR + Mclust 0.755 0.401 –
MNARz + SelvarMix 0.328 0.087 1 (missing variable)

The poor performance of the MNARz + SelvarMix baseline prompted a deeper investigation. To test whether properly handling imputation uncertainty could rescue the decoupled strategy, we also ran Multiple Imputation (MI) variants (missRanger-MI and MNARz-MI with random initializations) and pooled the results. Using the default Rmixmod backend for SelvarMix yielded the results below.

Table 7: Performance on Mixed (MAR+MNARy) Data (Rmixmod backend)
Method ARI NRMSE Relevant Variables
SelvarMNARz (Ours) 0.778 0.162 1, 2
Decoupled (SI MNARz) 0.328 0.263 1
Decoupled (MI missRanger) 0.336 0.257 1
Decoupled (MI MNARz) 0.328 0.331 1
Table 8: Performance on Mixed (MCAR+MAR+MNARy) Data (Rmixmod backend)
Method ARI NRMSE Relevant Variables
SelvarMNARz (Ours) 0.773 0.164 1, 2
Decoupled (SI MNARz) 0.328 0.387 1
Decoupled (MI MNARz) 0.336 0.322 1
Decoupled (MI missRanger) 0.328 0.248 1

From Table 7 and Table 8, we deduce two key takeaways: (i) Our joint model retains high ARI and correct selection under mixed mechanisms, including MNARy misspecification; (ii) MI does not repair the decoupled pipeline in this setting, with results mirroring single-imputation. This suggests two contributing factors:

  1. 1.

    Uncertainty Propagation. Decoupled pipelines treat completed data as observed, discarding posterior uncertainty in 𝒚m. Our EM-based framework propagates this uncertainty through all parameter and role updates.

  2. 2.

    Downstream Stability. The variable-selection backend (SelvarMix with its Rmixmod engine) shows instability (e.g., sensitivity to local maxima); MI then averages multiple weak fits. The joint estimation procedure is empirically more stable in this context.

Finally, under these mixed mechanisms, the class-level MNARz parameters in our model often converge to similar values across clusters for variables whose missingness is effectively MAR or MCAR, while remaining discriminative where missingness is truly class-linked. This helps explain the model’s robustness to misspecification.

Effect of Initialization on EM Stability and Accuracy. We compared a hierarchical clustering (HC) based initialization (Ward’s linkage on Euclidean distances with cluster centers extracted from the dendrogram cut at K) against multiple random initializations (MIs). HC places initial centers in high–density regions, yielding more stable responsibilities at the first E–step and fewer poor local optima than purely random starts. As shown in Table 9, HC attains higher ARI on 7 of 8 scenarios, and trails slightly once, confirming its overall robustness and improved convergence behavior for the EM algorithm.

Table 9: Comparison of EM initialization methods under 8 data scenarios in Section 5 (higher ARI is better). MIs: Multiple random initializations.
Scenario MIs (ARI) HC (ARI)
1 0.283 0.317
2 0.479 0.533
3 0.551 0.536
4 0.441 0.654
5 0.760 0.771
6 0.711 0.778
7 0.716 0.780
8 0.785 0.786

HC initialization offers a more sensible and stable warm start than multiple random restarts, typically improving both EM convergence and final clustering accuracy (ARI), with negligible overhead relative to the overall EM cost.

F.6 More on Transcriptome Dataset

Background and Prior Analyses

Dataset. We analyze the Arabidopsis thaliana transcriptome comprising 1267 genes measured across 27 experimental conditions aggregated from seven projects P1-P7. Genes were preselected for differential expression at least once in the hypocotyl growth switch time course (Project 6), making P6 biologically central. Following [37, 39], we retain all 1267 genes: 1149 are fully observed; 118 contain missing entries (107 with one, 10 with two, 1 with three). Overall, 9.3% of genes have any missingness and the global missing rate is 0.38%.

Prior findings. SelvarClust [37] (complete cases) found that including irrelevant variables degrades homogeneity; variable selection produced more coherent clusters and recovered known co-expression groups (e.g., a cluster of 15 genes co-clustered with 4 well-studied markers). SelvarClustMV [39] (MAR setting) expanded the gene set by reprocessing previously excluded genes and treating missing data within an EM framework, concluding that P6 (hypocotyl switch) and P7 (isoxaben treatment) are clustering‐relevant, together with P1-P4, whereas P5 (nematode infection) is not primarily grouping. Both studies support P2 (iron signaling) as a core axis for defining co-expression groups. Note that the 2012 analysis assumes MAR for the missing entries. In contrast, our framework explicitly models class-dependent missingness (MNARz), while retaining MAR/MCAR as limiting cases through parameterization.

Additional Interpretation of Our Results on the Transcriptome

Global outcome. Fitting SelvarMNARz for K∈2,…,20 with c=5, spectral distance weights 𝑷k, and pk​L​C structure, the selected model yields 18 clusters and a global role assignment with P1-P4 in 𝕊 and P5-P7 in 𝕌. This agrees with prior work on P5 (not clustering-relevant) but differs by reclassifying P6-P7 from 𝕊 (in [37, 39]) to 𝕌.

Cluster-level diagnostics clarify the difference. Let Rk2 denote the coefficient of determination from regressing the 𝕌-block (here, P5-P7) on the 𝕊-block (P1-P4) within cluster k. Table 2 shows pronounced heterogeneity:

  • •

    Several clusters (e.g., 6, 7, 8, 10, 12, 18) have Rk2>0.60, indicating that P5-P7 are largely explained by P1-P4 within those clusters. For these groups, assigning P5-P7 to 𝕌 is appropriate.

  • •

    A large aggregate (Cluster 1) exhibits low R2 alongside flat 𝕊-profiles and detectable 𝕌-activity, suggesting local decoupling of P5-P7 from P1-P4. This mirrors the biological intuition behind earlier inclusion of P6-P7 in 𝕊.

Hence, both the earlier global 𝕊-assignment (P6-P7 in 𝕊) and our global 𝕌-assignment (P6-P7 in 𝕌) are locally valid—but on different clusters. The discrepancy is explained by heterogeneity: the relationship between early axes (P1-P4) and late/stress projects (P5-P7) is cluster-specific.

Why does our global BIC prefer P5-P7 in 𝕌? Two factors are at play:

  1. 1.

    Global model selection. BIC aggregates fit-complexity tradeoffs across all clusters. Because many clusters exhibit high Rk2 (P5-P7 explained by P1-P4), the global criterion prefers a parsimonious 𝕊 (P1-P4) and assigns P5-P7 to 𝕌.

  2. 2.

    MNARz identifiability and shrinkage. Under MNARz, class-specific missingness parameters ρk​d absorb class-linked absence patterns. When a project effectively behaves as MAR/MCAR for many clusters, the fitted ρk​d across k becomes nearly homogeneous and the conditional dependence of P5-P7 on P1-P4 becomes tighter, further favoring 𝕌 globally.

Refer to caption
Figure 9: Mean expression profiles for each of the 18 clusters, displayed in separate panels. Light region indicates irrelevant P.

Biological reading consistent with both views. P2 (iron signaling) is reaffirmed as a core driver. P5-P7 behave as late/stress outputs strongly coupled to P1-P4 in many clusters (high Rk2), but retain independent variation in a sizable subset (e.g., Cluster 1). This reconciles the prior decision to place P6-P7 in 𝕊 (emphasizing hypocotyl-centric signals) with our global 𝕌 assignment (emphasizing aggregate parsimony across all clusters).

Practical implication. Our diagnostics suggest a natural extension: cluster-adaptive role assignment or a hierarchical prior tying per-cluster roles, which would allow P6-P7 to enter 𝕊 only where local evidence (low residual error) warrants it while retaining parsimony elsewhere. This aligns with the heterogeneity revealed by Rk2 and preserves the strengths of both global viewpoints.

Our unified 𝕊/𝕌 decision is globally efficient and statistically supported by MNARz-aware likelihood, while cluster-level diagnostics uncover biologically meaningful deviations. Together, they provide a coherent picture: P1-P4 are the principal axes; P5-P7 are predominantly redundant but locally informative in specific clusters, explaining the divergence from prior MAR-based analyses.

Appendix G Additional Details on Related Work in the Literature

G.1 Model-based Clustering

Model-based clustering conceptualizes clustering as a statistical inference problem, where the data are assumed to be generated from a finite mixture of probability distributions, each corresponding to a latent cluster. Unlike heuristic-based methods, this paradigm enables principled inference, allowing for parameter estimation via likelihood-based techniques and objective model selection to determine the number of clusters.

A prototypical instance is the GMM, wherein each cluster is characterized by a multivariate Gaussian distribution. Parameter estimation is typically conducted using the EM algorithm [16], yielding soft assignments in which each observation is associated with posterior probabilities across clusters. GMMs offer more flexibility than simpler methods like k-means, as they can model clusters with varying shapes and overlapping regions via covariance structures.

Bayesian mixture models (MMs) extend this framework by treating both the parameters and the number of components as random variables with prior distributions. This fully Bayesian approach enables comprehensive uncertainty quantification over both model parameters and cluster allocations. For a fixed number of components, inference is often performed using Markov Chain Monte Carlo (MCMC) or variational methods. However, Bayesian MMs face the label switching problem due to the symmetry of the likelihood with respect to component labels, which renders the posterior non-identifiable without additional constraints. Strategies such as relabeling algorithms and identifiability constraints have been proposed to address this issue [65]. A notable advantage of the Bayesian approach is the incorporation of priors on model complexity, e.g., via Reversible Jump MCMC or birth-death processes, to infer the number of clusters. Despite their flexibility and robustness, fully Bayesian clustering methods can be computationally demanding, particularly when MCMC chains converge slowly or require extensive post-processing to resolve label ambiguity [29].

G.2 Variable Selection for Model-based Clustering

Historically, variable selection was performed using best-subset or stepwise selection approaches, typically guided by information criteria such as AIC [1] or BIC [58]. While effective for moderate-dimensional settings, these methods become computationally prohibitive as the number of features D increases.

The advent of penalized likelihood methods improved scalability and enabled sparse modeling. The LASSO [70], a seminal technique using an L1 penalty, enables both variable selection and coefficient shrinkage. Predecessors include the nonnegative garrote [9] and ridge regression [28], the latter using an L2 penalty, which does not induce sparsity. Efficient algorithms such as LARS [18] and coordinate descent have facilitated high-dimensional applications of these methods. However, the time complexity of LARS limits its stability when D is very large.

Law et al. [32] proposed a wrapper approach that jointly performs clustering and variable selection through greedy subset evaluation, introducing the notion of feature saliency to assess variable importance in cluster discrimination. Andrews and McNicholas [3] developed the VSCC algorithm, combining filter and wrapper strategies: variables are ranked based on within-cluster variance from an initial clustering, then incrementally added subject to correlation thresholds to enhance separability. Their R package implementation, vscc, offers efficient noise filtering prior to model-based refinement.

Raftery and Dean [54] introduced a model selection framework for variable selection within GMMs, categorizing variables as clustering, candidate, or noise. Each candidate is evaluated using BIC to compare models where the variable does or does not influence clustering. Their framework accommodates statistical dependence between irrelevant and clustering variables via regression, avoiding overly simplistic independence assumptions. Scrucca and Raftery [59] enhanced this method, introducing computational heuristics in the clustvarsel R package to streamline EM evaluations.

Expanding this idea, Maugis et al. [38] proposed a three-role framework, relevant, irrelevant, and redundant variables, with redundancies modeled via linear regression on relevant variables. This design improves accuracy in scenarios with correlated predictors. Nevertheless, the stepwise search remains computationally demanding as dimensionality increases.

To address scalability, Celeux et al. [13] proposed a two-stage approach: first, variables are ranked by penalized likelihood using the method of Zhou et al. [75], which penalizes component means and precisions. Then, a linear scan through this ranking assigns roles. This heuristic drastically reduces computation, preserves model interpretability, and achieves strong empirical performance, though theoretical guarantees for recovering the correct 𝕊​ℝ​𝕌​𝕎 partition remain open.

Extensions to categorical data include adaptations of the Raftery-Dean method to Latent Class Analysis (LCA) by Dean and Raftery [15], with further refinement by Fop et al. [20], resulting in improved class separation in clinical datasets. Bontemps and Toussile [8] considered mixtures of multinomial distributions, using slope-heuristic-adjusted penalization to improve model selection under small-sample conditions.

G.3 Handling Missing Data

In practical datasets, missing values often occur and are categorized as MCAR, MAR, or MNAR. Under the MCAR assumption, missingness is unrelated to any data values; MAR allows dependence on observed data; MNAR, the most complex case, involves dependence on unobserved or latent variables.

Under MAR, Multiple Imputation (MI) [57] has emerged as a robust strategy to address uncertainty. Instead of single imputation, MI generates multiple completed datasets via draws from predictive distributions, followed by Rubin’s combination rules to aggregate inference. MICE [71], a flexible implementation, fits univariate models conditionally and iteratively imputes missing values, supporting mixed data types and nonlinear relationships.

High-dimensional settings pose new challenges, where full joint modeling becomes unstable. To address this, Zhao and Long [74] incorporated Lasso and Bayesian Lasso into the MICE framework to enhance prediction accuracy through regularized imputation, effectively reducing overfitting and variable selection bias in high-D regimes.

Machine learning techniques have also advanced imputation. MissForest [63] uses random forests to iteratively impute variables based on others, capturing nonlinearities in a nonparametric manner. Although not explicitly designed for clustering, the ensemble mechanism implicitly reflects local data structures akin to clustering. Deep learning methods, such as MIWAE [36], adopt autoencoder architectures to generate multiple imputations from learned latent spaces, enabling further extensions like Fed-MIWAE [4] for privacy-sensitive contexts. These models assume a latent manifold structure and are best suited for large sample sizes, though interpretability remains a challenge.

Addressing MNAR data in clustering models remains difficult due to identifiability issues and the need for strong assumptions or auxiliary information. Two general strategies exist: (1) Selection models, specifying a joint model for data and missingness mechanisms, and (2) Pattern mixture models, defining distributions conditioned on missingness patterns and modeling their influence on cluster membership. In [62], Sportisse et al. investigated identifiability conditions in MNAR selection models and proposed methods to augment clustering models with missingness-informed constraints to improve identifiability.

Cite this paper

Please cite the published version. Venue: NeurIPS 2025, Advances in Neural Information Processing Systems 38 (2025). DOI: NeurIPS 2025 Proceedings. Official record: NeurIPS 2025.

BibTeX
@inproceedings{ho2025unified,
  title     = {A Unified Framework for Variable Selection in Model-Based Clustering with Missing Not at Random},
  author    = {Ho, Binh H. and Nguyen Chi, Long and Nguyen, TrungTin and Nguyen, Binh T. and Hoang, Van Ha and Drovandi, Christopher},
  booktitle = {Advances in Neural Information Processing Systems 38 (NeurIPS 2025)},
  year      = {2025},
  url       = {https://openreview.net/forum?id=n7yVbKH7c3},
}