A Unified Framework for Variable Selection in Model-Based Clustering with Missing Not at Random
Authors: Binh Huu Ho
Long Nguyen Chi
TrungTin Nguyen
Van Ha Hoang
Binh T. Nguyen
Christopher Drovandi
Venue: NeurIPS 2025 poster ()
Keywords: Variable Selection
Finite Mixture Models
Model-Based Clustering
Missing Data
Asymptotic Selection Consistency
Primary area: general_machine_learning
TL;DR: Unified, provably consistent model-based clustering that jointly selects variables and handles MNAR via a data-driven penalty and explicit missingness–class modeling, validated on transcriptomics.

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.

A Unified Framework for Variable Selection in Model-Based Clustering with Missing Not at Random
Binh H. Ho*.1.2 Long Nguyen Chi+.3, Trung Tin Nguyen Souyen*.t.4.. *,1.4,5
Binh T. Nguyen Van Ha Hoang!.2, Christopher Drovandi Faculty of Mathematics and Computer Science, University of Science, Ho Chi Minh City, Vietnam.
2 Vietnam National University Ho Chi Minh City, Vietnam.
Schooo Information and Communications Technology,
Hanoi University of Science and Technology, Ha Noi, Vietnam.
ARC Centre of Excellence for the Mathematical Analysis of Cellular Systems.
Schoool Mathematical Sciences, Queensland University of Technology, Brisbane City, Australia.
Abstract
Model-based clustering integrated with variable selection is 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 preva- lent 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 unified framework designed to address these challenges simultaneously. Our approach incorporates data-driven penalty matrix into penalized clustering to enable more flexible variable selection, along with 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.
Introduction
Model-based Clustering. Model-based clustering formulates clustering probabilistic infer- ence task, assuming data are generated from a finite mixture model, with each component corre- sponding to cluster. This enables likelihood-basee Gaussian mixture models (GMMs) [42, 43] are a classical instance, which can be estimated via (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, 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
estimation and principled model selection.
the sopectation-maxiniin
"Co-first author, 'Coresponding author.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

label switching persist [64]. Bayesian nonparametric approaches, including reversible-jump MCMC and birth- 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 informa- tion criterion (AIC) [1] and Bayesian information criterion (BIC) [58], but become computationally infeasible as dimensionality increases. Penalized likelihood approaches, notably least absolute shrink- age 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 step- wise method for identifying clustering- relevann variables, later optimized by Scrueca [59]. Finally, Maugis et al. [38] enrich this framework by incorporating redundancy modeling, enhancing selection accuracy in correlated settings 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 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 single linear pass through the ranked list, replacing combinatorial search with faster and interpretable procedure that maintains competitive clustering accuracy. However, there is no theoretical guarantee that this ranking recovers the true 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 principled framework by replacing each missing value with several plausible alternatives, then combining analyses across imputed datasets account for uncertainty [57]. The multivariate imputation by chained equations algorithm [71] is a flexible implementation that fits sequence of conditional models. accommodating mixed data types. In nighdiiiensional 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-base 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-weightee autoencoder [36] and its federated variant [4), rely on latent variable representations to generate imputations, assuming data lie near a www-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 unobserved values. They showed that identifiability can be achieved under structural assumptions, such as requiring each variable's '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.
2

---

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 for recovering the true signal-redundan conduundant-uiinfommatvy partition. Our method achieves substantial empirical improvements in Section over existing approaches. This yields a unified, high-dimensional model-based clustering method that jointly addresses variable selection and MNAR inference in 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 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 summary, limitations, and future directions in Section Proofs and additional details are provided in the supplementary material.
Notation. Throughout the paper, we use the shorthand [NV] denote the index set 21. N for any positive integer N E N. We denote It as the set of real numbers. The notation [S] represents the cardinality of any set S. For any vector
RD, || denotes its p-norm. The operator O indicates the element-wise (Hadamard) product between two matrices. We use n E [N] to index observations and d E [D] to index variables. The vector yn E R denotes the full data vector for observation n, which may be further specified in more granular forms. The binary vector C2 = (CnI..... CDD {0, 13 represents the missingness mask, where Cnd = 1 if ynd is missing. The latent variable in E [K] denotes the mixture component assignment for observation n, with indicator variable 2nk: = 1 [2n=R)- We use E(-) 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 probabilistic framework by assuming that the data matrix Y = (y)..... yN) E R NDD with each yn E RR, is independently drawn from location-scale mixture distribution, 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 yn under a K-component GMM with parameters a a= [T, 14, Σ} given by: faMM K, a.) = Σk= TT. (gy 14g, Σ₁) where Ty > 0 and Σx= & = The term Nggn [A:: Σj.) denotes the multivariate normal density with mean vector N& and covariance matrix Σx. which is compactly denoted by O₂ = (1rR: Σg). The covariance matrix Σn 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)) 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-Whit 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 Model Selection Problem in Model-based Clustering, In high- -igeensimni 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 SRUW model, extending their earlier work [37], which assigns variables to one of four roles. Let denote the set of relevant clustering variables. Its complement, 3 comprises the irrelevant variables and is partitioned into subsets U and W. Variables in are linearly explained by a subset C S, while variables in W are assumed independent of all relevant variables. This framework facilitates variable-specine interpretation and avoids over- wver-enalizzaton from complex covariance structures, as typically encountered in constrained GMMs [12]. The joint density under the SRUW model, combining mixture components for clustering, regression, and independence, is defined as:
(aR.wwy) K,m r, 1, V, θ)= = Salaay K,m, a:) (y r, a y, 52) Simtep (y 11,7, D). (I)
3

---

Here, V = (S, U, W) denotes the variable partition, and @ is the full parameter set. The com- ponents are defined as: fels == fCMM. ( a 3, .2) == Ny a 3. Ω), with E R R 3 R [ and 52 E R|U structured by r. The independent part is fmdep = N, with variance structure / and covariance r. Model selection proceeds by maximizing BiC-type criterion:
(2) over (K, m. r. 1, V), where each term scores the corresponding model component, Although this ap- proach offers fine-grained variable treatment, it can be computationally demanding in high dimensions due to stepwise selection for clustering and regression.
critsic K, m, 1, clCalust K, + BIC (Y" r, Y²) + BICindep (Y (t),
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 Pk for adaptive regularization in Gaussian graphical mixture models. This method improves upon stepwise procedures [38, 37J by handling high-dimensional data more efficiently. The penalized log- ligeliiei is given by:
K
Σ In
A Σ
Σ
( ④ E )adl
k=1
k= diged
Here, Yn denotes centered and scaled observations, and Ps is typically chosen be inversely propor- tional to initial partial correlations, enabling adaptive shrinkage akin to the adaptive lasso [77]. Given grids of regularization parameters 9X and Gp. the EM algorithm from [75] is used to estimate the parameters ã(\,,) p) = (*,, Tan p),.).
p). En (A, Σx (X, p)). For each variable E [D]. clustering scores (x(d) are computed, and the relevant set $ is formed from those that improve BIC, while U captures non-informative variables. The final model is selected by maximizing critaic M,m, r, V.VVVK V(K,,.)) in Equation (2), where ((K,m)) = (S(,,) m)))) (K),) W ((K.))).
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 SRUW under MAR. key aspect is expressing the observed-data likelihood in a tractable form. Under the MAR assumption, the SRUW model admits global GMM representation, extending the result of [39] for the SR model. Given the density Equation (1) and Gaussian properties, the observed likelihood becomes
N
K
NV (va ikk,o Ancu))
f(Y K, m,,,,rr,y V θ
II
Σ
r==1
Vko0 %
O
where i² =
Ák,O
Here, Vk Yo and AR,00 Too correspond to
T 0o
the block means and covariances from the SR 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 SRUW 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 SRUW. To explicitly address MNAR data, 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 C') are retained (see Section 4). Let yn. (yf. y y denote the full data for observation n, partitioned by role (Section 2), and let c, be its binary missingness pattern. With 2n E [K] denoting the latent class, the complete-data density under 2nk 1 is
f(yn, Cn. 2nk = 1; ax, 4/a)) = Selust ( axx) Free Ml yh Oreg)
Finde( (mdep) fNINARZ(C, 2nk = 1; xb4),

---

where fMNARR(C,, 2nk = 1; = IT= ( (1 p(tb))) 1-end and p(zbk) is the class- specific missingness probability. The model parameters are Q (T1,...., TK, {ork ) = (0rg: Omdep) J: { J 1). The observed-data log-likelihood is:
N
K
{ Y. = Σ
log
TTfk ag ())c sha),
(5)
n=1
1=1
with the MNARz mechanism factored out due to independence. The complete-data likeliiil is:
K
Coomp! Y,,,
Σ
Σ
Znk log f ax, E)fc(Cni sbx)..
(6)
n=
To accommodate both MAR and MNAR patterns, we extend the model with a binary partition of variable indices: DAAR U DMNAR = [D], [DMAR| DM. [PMNAR| DMr. Variables in DMAR follow a MAR mechanism, and those DANNAR follow MNARZ:
fu (Yn: 0.k) Jetus(efr fro)) y Oreg) (v (mdep)) (0 mmdep);
MMR (Gr,M (M whi)) = II Pa(yni vla nnd (1 Pa(y vaa) '-cnd
dEDMAR
MANAR (CD, - II p(abx) (1 -1-ple (1-(((b)) (xx)))=cmu 2-cnd
de DANAR
The parameter blocks expand θ (T1.-..,, RK, fax) LbM, (p+))=))) Let yn (wh y,"), and define Cn.M (Cnd) Co.M' (Cud)de Then, since the MNARz mask is independent of yn given k, the observed-dat likelihood becomes:
(0:Y C)
Σ log Σ
T Se
MNNRR
(Co, fR (w or tbas)
fMAR (Ca.M yh wbm) dy The MNARz term factors of the integral, preserving the ignorability property for that block [62]. Only the MAR component requires integration or imputation.
where fř M (y%: aus xbx) : fa(yni (4))
Adaptive Weighting for Penalty of Precision Matrix. To enhance performance, we replace the inverse partial correlation scheme with spectral-based computation m in Equation (3). Starting from the initial precision matrix (0) we construct an unweighted, undirected graph ( k.ij | Tadj- Let (k)
((k) g()), where V = [D]. and include edge (i,j) in (() if
and D denote the adjacency and degree matrices, respectively. The symmetrically normalized Laplacian is then defined as L = (D()) (1) ( which offers scale-invariant measure of connectivity. We aim to shrink v toward diagonal target (i.e., an empty graph), with
Lurget The spectral distance between + (0) and this target DLs (I (0)) = [|spec(L ll,, and
the corresponding adaptive weights are given by PR.If DLS ( (0) +
Parameter Estimation. Following [75], we adopt similar parameter estimation for variable ranking and focus here on estimating SRUW under the MNARZ mechanism using the EM algorithm. Other estimation details are provided in the Supplement. Let each observation be partitioned as = (y% y/" and define & (Oreg. Quddep). Since fo independent of yn. the observed-datt log- likelihood becomes:
(O:Y" =
log
Tk Telnns bh))) on.. E),
n=1
1=
where (1 ag. ) = f felua (y
OKk Freg !
F. Grey Sintep
1 Oindep) dy
The E-step computes responsibilities:
(t-1) g(t- (111) [ (t-1) fo (h:: a
(t)
Σk मो
(t-1) (y:: att (1-1) 1) (cnn (1()-
5

---

The Q-function to be maximized in the M- M-step is:
Q(0 Σ
nk (log Th + g, (y(()) + go(bba)).
99 ( g flog fa (3) ar Ex) yn Cuz 2ndk = 1]
(e((bk) B 13 (log Se(Cn yni 1ba) yn Cr. Fick =
In the M-step, update:
N
nk
(1)
(1)
arg max
Σ
(0) nkgy (ax);
(akg)(ax)
arg max
(bbk)
5()) = arg max
(e) kkgy
n2k
We now summarize the workflow; the penalty acts only in Stage (ranking), while Stage B performs unpenalized SRUW model selection by single pass over the ranked list.
Algorithm 1 High level two-stage BRUW-MNARS procedure
Input:
incomplete data Y, mask C, model grid (K, m): regularization grids Gx (and Gp). Stage A: Ranking (penalized GMM fast-imputed data)
Produce fast single imputation y only to enable ranking.
For each (A. p) 9x x Gp. fit the penalized GMM to Y where Y denotes the centered/scaled Y: record p).
Compute ranking score Ox(d) for each variable d by counting along the path how often {Ara(A, p)] remain nonzero; sort variables by Ox (d) (see [13] for details).
Stage B: Role assignment (SRUW on original incomplete data)
Traverse ranked variables once (similar to [13]); at each step, fit unpenalized incomplete Y and decide the role (S/U/W) of the new variable using BIC criterion using Equation (2) as in [38]
Return ($, RR U, W), 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 robust warm start: (i) fast single imputation Y (only for ranking); (ii) K meeantt//heerrchhcall clustering/random start Y for ( 14x, Es); (iii) given initial partitions from step (ii) v set to diagonal estimates using samples assigned to each cluster (cluster-aware estimates); (iv) zha initialized from class-wise missing after soft E-step. We then
X-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 X, with the hard partition z from initialization.
Amax
max kE[K]]
(0) x)kd
A =
Amax
Aannn Amnaxe
Amin = Edmax: = [Z]]
This matches the KKT threshold which the l penalty zeros all μ4 entries. For let S₃. be the (hard-label) empirical covariance and note our Stage-A glasso objective uses the coefficient (2p/ms) in front of the weighted & term. A diagonal-forcing threshold is
T&k 1(5))
Pmax
max max &
$$1
(Px)
followed by the same T--point geometric path Pe: Pmaxx (Pmin/ Pmax) ) ((-1)/(L-1) with Pmin = EPpax We set diag P&) 0 when we do not impose penalty 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 SRUW model. We also show that these guarantees apply 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 (S, IR, U, W), let 0 (..R.U,))
denote the true parameter and 0 (,,) its maximum likelihood estimator.
Theorem (Informal: Identifiability of the SRUW Model). Let (K,,, and (K*, "",""""""""""""" ,"",""""""""""""""""""""" denote two models under the MNARZ mechanism. Let (K......)) S Y (K.....1.)) denote the parameter space such that each elemens O = (a, a, 3. 7, I) satisfies the following: (i) the component-spect
parameters (JuK, Σk) are distinct and satisfy. for all 8 S S, there exist 1 S k < k I. K such that pk,55s + *Mk, or Ek,a's i + Ek'als or Ek.aal + Ew salas where 3 denotes the complement of in S; (ii) I + 0. then for every j E R. there exists u E U such that Bug # O, and each E U affects at least one E E with Buz # O; and (ii) the parameters T and 7 exactly respect the structural forms and I, respectively. If there exist O E O (K.m.. and 0* E O "m""r" d' 1")3)) such that f(- K, r,,, V = =f( K",m,, r ,,, 0*), then (K, m. r, 1,, = (K*, m rr ,,,,,,,"
,,,,,,, V and e =
Theorem 2 (Informal: BIC consistency for the SRUW model). Assume the MNARz mechanism, the regularity conditions specified in the Supplementary Material, and the existence of unique tuple (Ko, mo, ro, lo, Va) such that the data- data-generating distribution satisfies = [(- 0 ma No)) for some true parameter 0 where h is the density function the sample Y and the model specification (Ko, mo, To, la Vo) is assumed to be known. Let (Ko, mo, ro, la) be fixed. Then the variable selection procedure that selects the subsets (S, B. Ü W) by maximizing the BIC criterion is consistent, in the sense that IP
(S, RR Ú, W) (So. Ro, Wo, Wo))
N
Theorem (Informal: Selection consistency of the two-stage procedure). Under certain regularity conditions in Supplementary Material, the backward stepwise selection procedure in the SRUW model, guided by the variable ranking stage. consistently recovers the true relevant variable set S
with high probability. That is, PPS =Si) 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 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 L (ax) denote the empirical objective function based on sample size N, and let a' be the true parameter. There exist constants 71, ?2 > 0 such that for every perturbation in a restricted neighborhood of ax, (x(a** + A) [x((a*)= (V A) Z ( A) is tolerance term that decays at rate N 1/2, complete and rigorous proof is provided in the Supplementary Material.
- 1 + 7 | l TN A),, where TN
Implications of the theorems. In the framework of [13], the parameter e in the BIC-based selection step (representing the number of consecutive non-positive BIC differences allowed before terminating inclusion into S or W) 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 lD) (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 Ox (d) scores between the last true signal and the first noise variable is sufficiently pronounced, the exact value of e 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 e could lead to the inclusion of noise variables, especially if the variable ranking is imperfect.
7

---

5 Experiments
Our primary objective in this section is evaluate the clustering quality, imputation aceuracy, 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-impuuationnccussrring pipeline with the imputation method to be Random Forest [10] using missRanger package [40, 63J. 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 in [37, 39, 13] and interpret the result in comparison with previous findings in [37, 391. 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 for variable/model selection with missing values treatment benchmarking. The first data set mimics [39] where each observation yn E R is drawn from four-component Gaussian mixture on the clustering block $ = (yn, y2, ys). The components are equiprobable and have means 00,0.0) (-6, 6, 0), (o, 0, 6) (-6, 6,6) with a diagonal covariance matrix Σ = (aag(6v2,1,22 oaaak The redundant block U = {ys. ys is generated by a linear regression U = (-1,2) + Xs((0.5, 2) (1, + c, where R(0) denotes the usual planar rotation. The remaining W variables are iid N/O, 1) noise. In the second design, an observation yn E R The clustering block now consists of (y1 y2). drawn from equiprobable four-component Gaussian mixture with means (0, 0), (4. 0), (0. 2), (4.2) and common covariance 1. Conditional on (31, y2). the vector 3/3:14 follows a linear model 33:14 a+ (y y2) where a. the x 12 coefficient matrix B and the diagonal covariance 52 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 912:14- In our setting, we consider K E {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 2 and compute PA with the spectral distance. Further details of the synthetic-data generation scheme are provided in the supplement. From Figure we observe that SelvarMNARZ delivers the highest ARI and the lowest WNRMSE across every missing-date level in both scenarios. Its performance declines only modestly as missingness increases, while the theuultte 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 significantly slower performance decline under high missingness, we conducted one-sided Welch's 's t-tests (Bonferroni- correcctdddr 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.
2000
(1, 3) (TT(6))))) and
Table 1: Welch's 's t-test for ARI (MNAR, 50% missingness, a = 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
s * ** 8 * **
Selvar
0.348
0.108
vs. SelvarMNARz
< 0.001
VarSelLCM 0.344
0.101
SelvarMNARz
0.001
#*
While all methods show a monotonic performance drop missingness increases, these results confirm that the superior performance SelvarMNARz is not 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-date mechanism rate. On Dataset 2, every method succeeds un-
8

---

MAR
MNAR
MAR
MNAR
Dakeet
Dukiet_2
Dutaee22
Dusise_
Detaset22
Datset
DutAeeL2
0.00
WARNINE
Cuntonee
Himel
Clntonne
Seler
detar Hobestos
Clstaman
Castene
Figure I: Comparison of four models under MAR and MNAR mechanisms over 20 replications; for ARI/WNRMSE, higher/lower boxplots indicate better performance.
MAR
Cluster sumbers
Relewant vare
C Clurter sumbers
Relevant vars
024
188
1.2
0.5
MNAR
Dutaset rumbers
Dataset
2utaset.2
Dusiset vers
0.3
0.3
3.2
Model
Clsoaascl
Selver
Schu/MMARz
VasslLCM
Figure Proportions choosing correct relevant variables and cluster components over 20 replications.
der 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-gen geadeeeeeeenngggg designs: Dataset includes pure-noise variables W within the SRUW framework. whereas Dataset 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.
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 the Supplementary Material: here, we summarize the settings and key findings. We fitted SelvarMNARz for K (2..... 20} clusters with 5, spectral distance weights Pk, and LC mixture structure as [39]. For conciseness, we refer to Project as Pl, Project 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 broader gene set (including entries with missing
Transcriptome Data.

---

Dafcrere -F Prgeet 1-4, lgferose Projec 5-7
Figure Mean expression profiles across clusters. Light region indicates irrelevant
Table 2: Size of each mixture component and R obtained by regressing the U--bock the S-block within that cluster. Values extremely close to (marked by ") occur because the U--llok nearly constant in those groups.
Cluster # Genes
R²
Cluster
# Genes
R²
Cluster
# Genes
R²
715
0.02
7
36
0.68
13
46
0.57
11
1.00*
8
13
1.00"
14
19
1.00"
1.00*
15
1.00"
15
1.00*
0.45
29
0.59
16
1.00*
5
124
0.12
11
71
0.43
17
24
0.90
1.00*
12
24
0.91
18
10
1.00"
data) considered. Additionally, three methods agree that P5 is not primary clustering factor and is more likely a correlated effect explained by broader cellular processes. However, key departure lies the reclassification of P6 and P7, previously identified as informative clustering drivers [37, 39], into redundant U block, explained by P1-P4. This suggests that, after accounting for MNAR effects, variation P6-P7 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 R² 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 Rt2 > 0.60, indicating that late-stress variation is well encoded by P1-P4. Cluster 1, although large, is transcriptionally flat across both S and U; its low R 0.02 thus reflects vanishing denominator rather than hidden structure. True residual structure (unexplained U signal) is concentrated 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- class dependencies. This enhances clustering flcxibility and robustness in nigh-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 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 statistically grounded stepwise approach. Bontemps and Toussile [8] proposed mixture of multivariate multinomial distributions with combinatorial variable selection and slope heuristics for penalty calibration. Building upon these, 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.
10

---

Acknowledgments
TiungTin 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 express their special thanks to the authors of [37, 39], particular Professor Cathy Maugis- Rabusseau and Dr. Marie-Laus Martin, for providing the Arabidopsis thaliana transcriptome dataset.
References
[1] H. Akaike. new look the statistical model identification. JEEE transactions on automatic control, 19(6)711-723, 2003. 2, 70
[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-3) -132, 2009. 30
[3] J. L. Andrews and P. D. McNicholas. Variable selection for clustering and classification. Journal of Classification, 31(2): 133-153, 2014.
[4] Balelli, Sportisse, Cremonesi, P.-A. Mattei, and M. Lorenzi. Fed-miwae: Federated imputation of incomplete data via deep generative models. arXiv preprint vrXiv:2304.08055.
2023. 2,
[5] E. Barut, J. Fan, and A. Verhasselt. Conditional sure independence screening. Journal of the American Statistical Association, 1 (2651)):126612215 2016, 2
[6] J.-P. Baudry, C. Maugis, and B. Michel. Slope heuristics: overview and implementation. Statistics and Computing, 22(2):455 470. 2012.
[7] C. Biernacki, G. Celeux, and G. Govaert. Assessing a mixture model for clustering with integrated completed likelihood. IEEE Trunsactions on Puttern Analysis and Machine Intelligence, 22(7):719-725, 2000.
[8] D. Bontemps and W. Toussile. Clustering and variable selection for categorical multivariate data. Electronic Journal of Statistics, 27234-23771 Jan. 2013. 10, 71
[9] L. Breiman. Better subset regression using the nonnegative garrote. Technometrics, 37(4)3373- 384, 1995. 71
[10] Breiman. Random forests. Machine Learning, 45:5-32, 2001. 8
[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. 4, 23
[12] G. Celeux and G. Govaert. Gaussian parsimonious clustering models. Pattern Recognition, 28(5):781-793, 1995. Publisher: Elsevier.
[13] G. Celeux, C. Maugis- Rabusseaa,, 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. 2, 3, 4, 6, 7, 8, 24, 60, 61, 67,
[14] M. C. Chong. H. D. Nguyen, and TrungTin Nguyen. Risk Bounds for Mixture Density Estimation on Compact Domains via the h-Lifted Bellback--eibble Divergence. Transactions on Machine Learning Research, 2024.
[15] N. Dean and A. E. Raftery. Latent class analysis variable selection. Annals of the Institute of Statistical Mathematics, 02:11-35, 2010. 10, 71
[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).133,, 1977. 1, 70

---

[17] D. Do. L. Do, S. A. McKinley, J. Terhorst, and X. Nguyen. Dendrogram of mixing mea- sures: Learning latent hierarchy and model selection for finite mixture models. arXiv preprint arXiv:2403.0/
2024. 3
[18] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression. Annals of Statistics, pages 407-451, 2004. 2, 71
[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. 2
[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-21 10, 2017. 10, 71
[21] Forbes, 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. 3
[22] Forbes, H. D. Nguyen, T. Nguyen, and Arbel. Summary statistics and discrepancy measures for approximate Bayesian computation via surrogate posteriors. Statistics and Computing, 32(5):85, Oct. 2022.
[23] R. Foygel and M. Drton. Extended Bayesian Information Criteria for Gaussian Graphical Models. In Lafferty, C. Williams, Taylor,, Zemel, and A. Culotta, editors. Advances in Neural Information Processing Systems, volume 23. Curran Associates, Inc., 2010. 3
[24] Schhnatt--Sc C. Pamminger, A. Weber, and R. Winter-Ebmer. Labor market entry and earnings dynamics: Bayesian inference using mixtures-of-expert Markov chain clustering. Journal of Applied Econometrics, 26-))11 2012.
[25] C. R. Genovese and L. Wasserman. Rates convergence for the Gaussian mixture sieve. The Annals of Statistics, 28(4): 105 1127, 2000. 3
[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. 3
[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-30) 2016.
[28] A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics. 12(1):55-67, 1970. 2, 71
[29] A. Jasra, 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 -1::0066,, 2005. 2, 70
[30] S. Kim, M. G. Tadesse, and M. Vannucci. Variable selection in clustering via dirichlet process mixture models. Biometrika, 93(4):877-88) 893, 2006. 2
[31] J. Kwon and C. Caramanis. The em algorithm gives sample-optimality for learning mixtures of well-separated gaussians. Conference on Learning Theory, pages 2425-2487, PMLR, 2020.
62
[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-
1 1166, 2004. 2, 71
[33] R. J. Little and D. B. Rubin. Statistical analysis with missing data. John Wiley & Sons, 2019. 23
[34] P.-.. 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. 23
12

---

[35] M. Marbac and M. Sedki. Variable selection for mixed data clustering: A model- based approach. arXiv preprint arXiv: 1733002003 2017. 8, 23
[36] P.-. Mattei and J. Frellsen. MIWAE: Deep Generative Modelling and Imputation of Incomplete Data Sets. K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages
PMLR, June 2019. 2, 71
[37] C. Maugis, G. Celeux, and M.-I.. Martin Mannieete Variable selection for clustering with gaussian mixture models. Biometrics, 65(3):001-70, 2009. 3, 4, 8, 9, 10, 11, 24, 29, 68, 69
[38] C. Maugis, G. Celeux, and M.-L. Martin-Magniette. Variable selection in model-based cluster- ing: A general variable role modeling. Computational Statistics & Data Analysis, 53(1 1):3822 3882, 2009. 2, 3, 4, 6, 7, 23, 24, 27, 29, 40, 71
[39] C. Maugis-Ratie Rabusseau, M.-L. Martin- Martti-MMagnictte, and S. Pelletier. Sclvarclustmv: Variable selection approach in model- clustering allowing for missing values. Journal de la Société Française de Statistique, 153(2)221-36, 2012. 3, 4, 8, 9, 10, 11, 68, 69
[40] M. Mayer. missRanger: Fast Imputation of Missing Values, 2025. R package version 2.6.2, https://mayer:: watpaaerrrrggtttttu github..hu.... ii/iii/iiii io/misssanger//.......
[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. 62
[42] G. J. McLachlan and K. E. Basford. Mixture Models: Inference and Applications Clustering, volume 38. Marcel Dekker, 1988. 1
[43] G. J. McLachlan and D. Peel. Finite mixture models. John Wiley & Sons, 2004.
[44] S. Negahban and M. Wainwright. Restricted strong convexity and weighted matrix completion: Optimal bounds with noise. Journal of Machine Learning Research, 13(53):1665-1697 2012. 23, 47
[45] D. N. Nguyen and Z. Li. Joint learning of Gaussian graphical models in heterogeneous dependencies of high-dimensional transcriptomic data. The 16th Asian Conference on Machine Learning (Conference Track), 2024.
[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. 3
[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. 3
[48] H. D. Nguyen, Nguyen, and F. Forbes. Bayesian Likelihood Free Inference using Mixtures of Experts. In 2024 International Joint Conference on Neural Networks (IJCNN), pages 1 1-8, 2024.
[49] T. Nguyen, F. Chamroukhi, H. D. Nguyen, and G. J. McLachlan. Approximation of probability density functions via docation-scale finite mixtures in Lebesgue spaces. Communications in Statistics Theory and Methods,
5059, 2023. 3
[50] T. Nguyen, D. 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.
[51] T. Nguyen, H. D. Nguyen, F. Chamroukhi, and F. Forbes. A non-asymptotic approach for model selection via penalization high-dimensional mixture experts models. Electronic Journal of Statistics, 16(2):4742 4822, 2022. 3
13

---

[52] T. Nguyen, H. D. Nguyen, F. Chamroukhi, and G. J. McLachlan. Approximation by finite mixtures continuous density functions that vanish infinity, Cogent Mathematics & Statistics, (())17700611, Jan. 2020. 3
[53] X. Nguyen. Convergence latent mixing measures finite and infinite mixture models. The Annals of Statistics, 41(1):370-4000 2013.
[54] A. E. Raftery and N. Dean. Variable selection for model-based clustering. Journal of the American Statistical Association, 101(473):168-178, 2006. 2, 71
[55] A. Rakhlin, D. Panchenko, and S. Mukherjee. Risk bounds for mixture density estimation. ESAIM: PS, 9:220-229, 2005. 3
[56] D. B. Rubin. Inference and missing data. Biometrika, 63(3):581-592 1976. 23
[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. 2,
[58] G. Schwarz. Estimating the dimension of a model. The annals of statistics, pages 461-464, 1978. 3, 70
[59] L. Scrueca 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. 2, 8, 71
[60] W. Shen, S. T. Tokdar, and S. Ghosal. Adaptive Bayesian multivariate density estimation with Dirichlet mixtures. Biometrika, 100(3):623-44,, 2013.
[61] C.-Y. Sin and H. White. Information criteria for selecting possibly misspecified parametric models. Journal of Econometrics, 71(1):207-225, 1996. 3
[62] A. Sportisse, M. Marbac, F. Laporte, G. Celeux, C. Boyer, J. Josse, and C. Biernacki. Model- based clustering with missing not random data. Statistics and Computing, 34(4):135, 2024. 2,
4, 5, 23, 60, 67, 72
[63] D. J. Stekhoven and P. Bühlmann. missffress---o--prraaeee missing value imputation for mixed-type data. Bioinformatics, 28(1):112-1138 2012. 2, 8, 71
[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. 2
[65] M. Stephens. Dealing with label switching in mixture models. Journal of the Royal Statistical Society: Series (Statistical Methodology). 62(4):795-80)) 2000. 70
[66] N. Städler, P. Bühlmann, and S. van de Geer. l l1-penalization for mixture regression models. TEST, 19(2):209-255 256, Aug. 2010. 23
[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)6622617,, 2005.
2
[68] H. Teicher. Identifiability of finite mixtures. The Annals of Mathematical Statistics, 3264):1265- 1269, December 1963. 30
[69] T. Thai, T. Nguyen, D. Do, N. Ho, and C. Drovandi. Model Selection for coussian-gater Gaussian Mixture of Experts Using Dendrograms of Mixing Measures. arXiv preprint arXiv:2505.1305: 2025.
[70] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1)2267-288 1996. 2, 23, 71
[71] S. Van Buuren and K. Groothhis-Oudshoonn mice: Multivariate imputation by chained equations in Journal of statistical software, 45:1-67, 2011. 2, 71
14

---

[72] J. Westerhout, T. Nguyen, X. Guo, and H. D. Nguyen. On the Asymptotic Distribution of the Minimum Empirical Risk. Forty-first International Conference on Machine Learning, 2024.
3
[73] D. M. Witten, Friedman, and Simon. New insights and faster computations for the graphical lasso. Journal of Computational and Graphical Statistics, 20(4):892-900, 2011. 62
[74] Y. Zhao and Q. Long. Multiple imputation the presence high-dimensional data. Statistical Methods in Medical Research, 25(5):2021-208) 2016. 2, 71
[75] H. Zhou, W. Pan, and X. Shen. Penalized model-based clustering with unconstrained covariance matrices. Electronic Journal of Statistics, 3(none): 1496, 2009. Publisher: Institute of Mathematical Statistics and Bernoulli Society. 2, 4, 5, 23, 71
[76] M. Zhou et al. Global convergence gradient em for aver-parameterized gaussian mixtures. arXiv preprint 1xXiv:25066065844 2025. 62
[77] H. Zou. The adaptive lasso and its oracle properties. Journal of the American Statistical Association, 101(476):4418-422) 2006. 4
NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
We provided technical proof and experimental evidence to support our claims.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer:
[Yes]
Justification:
We included it at the end of the paper.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model sell-specification asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g. if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
15

---

The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't It acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proof's
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: The paper provided the full set of assumptions and complete proof in appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We provided experiment descriptions in the main paper and appendix. Further- more, we will release code implementation if the paper is accepted.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access to a hosted model (c.g., in the case
16

---

of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (c.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results. 5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We will publish the code when the paper is accepted or through an anonymous link if some Reviewers ask for it.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: //nips.ccm. ttapps.cc/pbllcc//eccccc gwwde///e/////eee/
onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provided details for experiments in the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
17

---

The experimental setting should be presented in the core of the paper to level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:
[Yes]
Justification: We conduct experiments with appropriately reported error bars, defined in statistically sound and reproducible manner.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma eтror bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (c.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
This information is included in our results.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurIPS Code of Ethics ttttpp:///euriii http://uurrppsccppuuuiitse wttp:cccc/ttttcciieeeeeeec ctt.eiiieeecccccccieeeeeeees
Answer: [Yes]
Justification:
Our paper follows the NeurlPS Code of Ethics.
18

---

Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws or regulations in their jurisdiction)
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Our work contributes to reducing the carbon footprint associated with training statistical and machine learning models, thereby promoting more sustainable computational practices. We do not anticipate any negative societal impact.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question: Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Answer: [NA]
Justification: Our work does not pose such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
19

---

12. Licenses for existing assets
Question:
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We have properly cited papers and resources used in our experiment.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: Our paper does not release new assets.
Guidelines:
The answer NA means that the paper does release assets.
Researchers should communicate the details the dataset/code/model as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
We do not have experiments involving crowdsourcing or research with human subjects.
Justification:
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
20

---

15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
We do involve crowdsourcing or research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
Declaration LLM usage
Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing. editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: LLMs have been employed solely for writing, editing, and formatting purposes. Their use does not affect the core methodology, scientific rigor, originality of the research.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy (https: Ceeeurips.cc/Coofeeeeeeeee for what should or should not be described.
21

---

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 Section A, with the aim of enhancing the reader's understanding of our key theoretical and methodological developments. Next, summarize the definitions and theoretical results for the SR and SRUW models in Section C. We then present the regularity assumptions and technical proofs for the main theoretical results in Sections C and D, respectively. Finally, we include additional experimental results and further discussion of related work in Sections F and G.
Contents
Introduction
2 Background
3 Our Proposal: Variable Selection in Model-Based Clustering with MNAR
4 Theoretical Properties
5 Experiments
Conclusion, Limitations and Perspectives
10
A More Details Main Challenges and Contributions
23
B Preliminary of SR and SRUW Models
24
Regular Assumptions for Main Thcoretical Results
28
D Proof of Main Results
28
D.1 Proof of Theorem
29
D.2 Proof of Theorem 2
31
D.3 Proof of Theorem
40
Detailed EM Algorithms for the Two-Stage Procedure
57
Additional Experiments
60
E1
Metric Details
60
F2
Assessing Integrated MNARz Handling
60
F3
Sensitivity on the choice
60
F4
Computational Times and Scalability
62
F55
More on Simulated Dataset
67
F.6 More on Transcriptome Dataset
68
G Additional Details Related Work in the Literature
70
22

---

G.1 Model-based Clustering
70
G.2 Variable Selection for Model-hased Clustering
70
G.3 Handling Missing Data
A More Details on Main Challenges and Contributions
Model-based clustering with variable selection, particularly using the sophisticated SRUW framework [38], offers powerful approach to understanding data structure by not only grouping data but also by assigning distinet 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 frontier. Specifically, the MNARZ mechanism (where missingness depends only latent class membership) offers 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 nigh-dimensional regression (e.g.. LASSO [70, 75, 11J), 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- likelihoooo 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, attempt to make several contributions to address the aforemen- tioned challenges, particularly focusing on providing a theoretical foundation for a two-stage variable selection procedure within the SRUW-MNAR? framework:
1. Though the extension of the reinterpretation result from [62] is straightforward, this signifi- cantly simplifies the analysis and provides way handle identifiability of SRUW under MNAR missingness, which potentially facilitates algorithmic design.
2. Building on the MAR reinterpretation and existing results for SRUW data model identifia- bility and MNARz mechanism identifiability, we establish conditions under which the full SRUW-MNARY model parameters (including SRUW structure, GMM parameters, regres- sion coefficients, and MNARz probabilities) are identifiable from the observed data. This important since identifiability is prerequisite for consistent parameter estimation.
3. To underpin the variable ranking stage, which uses LASSO- penalized GMM, we provide rigorous analysis of the penalized GMM estimator based on [66, 44, 34]: First, we establish non-asympiotic sup-norm bound on the gradient of the GMM negative log- likelihood at the true parameters, which 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 specific cone, which is essential for sparse recovery. These settings provide the necessary
23

---

high-ding digh-dimensional statistical guarantees for the penalized GMM used in the ranking stage, ensuring that its estimated sparse means are reliable.
4. Leveraging the parameter consistency of the penalized GMM, prove that the variable ranking score Ox (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. For the main theorem 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 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- ranking, especially in the context of MNAR pattern. It bridges the theoretical gap between the computationally intensive and a more practical nigh-dimennione strategy.
clustering (SRUW) that incorporates an initial LASSO-II
6.
While full optimality of is not derived, we provide a condition for choosing e to achieve target false positive rate in relevant variable selection, given the single-decision error probability. This moves beyond purely heuristic choices for in [13], paving the way for more adaptive or theoretically grounded choices of in practice. This offers a more principled approach to selecting this single hyperparameter in the algorithm, enhancing the reliability of the selection procedure.
B Preliminary of SR and SRUW Models
This part summarizes definitions and theoretical results of the SR and SRUW models in [37, 38). We start with the following two definitions:
Definition 1 (SR Model for Complete Data), Let Yn E RR be the n-th observation,for for 1,...,0 N. A model MgR (K, m, S, R) is defined by:
1. K: The number of mixture components (clusters).
2. m: The covariance structure of the GMM for variables in S.
3. S C D): The non-empty set of relevant clustering variables.
4. IR C S: The subset of relevant variables in used to explain the irrelevant variables
ge (1.....) D} S via linear regression.
The parameters of the model are = (ax, 3, 32),, where a are the GMM parameters for variables S. is the intercept vector, 3 is the matrix of regression coefficients for explaining o from y and Q is the covariance matrix of the regression component. The complete data likelihood for one observation yn is:
flyn; m, S, 0) = (y;; K, m, Or) (3. +. 3,
where S denotes the sub-vector of yn corresponding variables in S, and similarly for y and y felust is a K-component GMM density, and fwee is multivariate Gaussian density.
Importantly, SR model can be equivalently written as single K-component GMM for the full data vector Un E R [37]:
K
f(y:: K. m. R,
Σ T D(y Uk. AE)
(8)
1=1
where (-: v. A) is the multivariate Gaussian density with mean and covariance A. The parame-
ters (Vr, 5 are constructed from the original SR parameters (ar, a. 00.3,,5 52) as follows: Let A [S) x |Se matrix derived from 3 (which is [R] x (8a1). For a variable j E S and E St, Agt Byn if the j-th variable of S is the jr -th variable of RR and Ast = if the j-th variable of R.. Then, for each component k:
24

---

The mean vector Da R has elements:
pkij (a + A),
if variable j S if variable j S
Pkj
where Pg is the mean of the k-th component for variables in S.
The covariance matrix 4 DxD has blocks:
ERsi (EgA)
WjE,les ifj E E, i 1 fj
A
(ATEa))
(Ω+ EEA))
Definition (SRUW Model for Complete Data), Let yn E RD be the n-th observation. An SRUW
model MsRUW (K,m, m, r, L, V) defined by:
1. K: The number of mixture components.
2. m. The covariance structure of the GMM for variables in S.
3. r: The form of the covariance matrix Q for the regression of U on R.
4. l: The form of the covariance matrix I for the independent variables W.
5. V (S, R, U, W): A partition of the variables, where
S: Non-empty set of relevant clustering variables.
R C S: Subset of S regressing variables in U. (R -R= -0
0).
Uj Set of irrelevant variables explained by R.
W: Set of irrelevant variables independent of all variables in S.
SUtU W = (1..... D} and these sets are disjoint.
The parameters are 0 = (a, a, 3,,, 7, I), where a is the GMM parameters. The complete data density for one observation ya
f(yn K, m, r, V. (y: K, m, freg + yh findep (V y, r)
This can be written K-component GMM for the full data vector yn E R.:
K
f(yn: K, F.V,
( Ux Ax)
where Un and Â constructed from O. Let (v (S, UJ) part derived as the SIR model (where plays the role of (ig) Then:
A (8.1)) be the mean and covariance for the (3.1)
(1)
(8.3)) -7))
U/ =
(3,U)
Å
due to the independence of W from and UV (given the cluster k, which implicit in the construction of V (8,U)
using )
Similar to global GMM representation of SR model under MAR, we can obtain such global form of SRUW model under the same missingness pattern, which is illustrated by the following claim.
Proposition 1. Under MAR, SRUW has a natural global GMM representation with appropriate parameters.
25

---

Proof of Proposition 1. Let Ur = (y yn yn be the complete data vector for observation n E {I.... [L......) N. partitioned according to the SRUW model structure V (S, R, UJ W). The parameters are 0 = (ox, a, [3 7, I',, where a (T ... Wh: Pi...... MK .,....; Σx)-
The complete-data density for single observation yn given its membership cluster k: (denoted by 1), is:
fyn Zik = 1, V, 0) = Sclua pux Σx) (we + 3,3,S) Qt Simdep Un 7, T) The marginal complete-dan density for bn summing over latent cluster memberships is:
fyn V, 0)
Σ Th [y pix En) frog (v v. S2) findopyy ; r)]
(10)
Due to the independence of yn from (y yh)) given the cluster and its parameters 7, T are not k-specific, we can rewrite the term inside sum as follows: suppose that f& [x Su4.,//. y = & (yy: n p&g Σ&) -by' a + 3n St). This is the density of (y v.) given cluster k: Then, Equa- tion (10) becomes:
fly V, 0)
Σ Tt.ee
Sindee (Yn 3 T)
The term (yy 7, 7, I) factors out of the sum over because 7 and are not indexed by k. Let = (y., v²") be the partition of yn into observed and missing parts for observation The WWMMa) ), where M,, denotes the set of missing variable missing parts are y = (y SnM, Yn unbba yn
indices for observation n. The observed-data density for observation is obtained by integrating fyn; V, 0) over y::
f(U V.
=
/ [(yh Un ; V) dy.
K
fk sno, , SnM, Unan LrMM))
Yn
X Siddep
(ya,
who WnM,,,
yn
r) ;
dyn SnM,
dyn dy
UnM., WnM,
where On denotes the set of observed variable indices for observation We can swap the sum and the integral:
K
fly:
[ y. yn y. UnM...
Sno.
8nM.,
Uno..
X S Sideep
(yn
w0o. WrM,, Wnnm; y, r) dy
Un
SnM, dy UnM. WnM, dy,
Since f-(-) only involves variables in S and U, and fehst (577, I) only involves variables in W, and these sets are disjoint, the integration can be separated:
Σg Tk / falueno
y
SoM,
yn
unc., yn Untba) dy dy
UnM.
If
Sindep((3 (y
Wno, yn
WnM dy WOM,
Let fu(yn (3.u)na.)
[a 0n y snic. yn uno.. yn UnM. d SnM. dy UnM.. This is the marginal density of the observed parts of (S, U) variables for component k. It is Gaussian, (8.U)m0 (3) (8.)) where these parameters are derived from the complete-dar data SR hoo
k,00
parameters for the (S, U) part. Let (yy x( r(2)) = Tyy wno, y WOM., ;7, I') dy W-M.,
This is the marginal density of the observed parts of W variables. Then,
fly:: V. 0)
Σ f (1)
(B.U)nD (8,U)(m)
((uu)())) R,00
finoe (v w00., ://00 (n) I!m)) (11)
k.o
26

---

From Equation (11), we can get
N
K
II
Σ
Th fchust D A (m)
koo
k,00)
n=l
by allowing Syn (y :0o 76")
r( to "absorbed" into the sum over So, we can write:
fly:: v,
Σ Th F ( (3,)0n (8.U)(n)
().)))))) k,o 2,00
finde((y) wro 10:" )
(())
Since the variables (S. U) and W are conditionally independent given k, the product of their marginal observed densities is the marginal observed density of their union. Let y
4 (8.
(yn
(,,
Let ((()) yn yn
WnD.)..
(1 (r)"))r)² Let
(n)
K000
koo
2,00
r (n)
Then, felust yn
(S,U)nO,.
(s,))m)) 8,00
(.))))))
finopp yn R0n 76"
(n)
I( felust (n koo
(n))
hoo
Therefore, the observed-data density for single observation n. is:
K
flyh V 0)
Σ
Th Telus Sh (m) koo A Mooo
(12)
And for N i..d observations, the total observed-data likelihood is:
K
II K, m=
II
Σ
Setus (1h: (1 (ndee)
(13)
When the parameters meet some distinct properties, SRUW under complete-data achieve identifi-
ability, as shown in Theorem 1.
Fact 1 (Theorem in [38]). Let (K.m,rd.)) be subset of the parameter set Y (KKm,.r.)) such that elements = (ax a, B, 22 y, r)
contain distinct couples (Mk, Σx) fulfilling Vs S ((k,
W K;
Pk,5|s + Mkl or Dx,5 + EE Bis Ek, # En ,51|
(14)
where denotes the complement in S of any nonempty subset of
U + 0.
for all variables j of R, there exists a variable of U such that the restriction Bug of the regression coefficient matrix associated to j and a is not equal to zero.
for all variables u of U, there exists variable j of B. such that Buj + 0.
parameters Ω and exactly respect the forms and I respectively. They both diagonal matrices with at least two different eigenvalues if r LB) and = [LB] and QL has at least non-zero entry outside the main diagonal ifr [LC].
Let (K, m,r,l, V) and (K* m², 1 W*) be models. If there exist O ① (KKm..)..)) and 0² E (ggr mm xr (*)3-) such that
[(JK, m. r, VV - /(-K"" m", ,,",r",1" 0")
then (K,m, m, 1,,, (,V) (K", m² r" r",,",",",""" v*) and 0 = 0" (up a permutation of mixture components).
27

---

C Regular Assumptions for Main Theoretical Results
For fixed tuple of (Kg, mg, To, lo: Vo), we write f(y; instead of fl; Ko, ma, So Ro. 0) or fy: V, for short notation.
Assumption (SR). There exists unique (Ko, mo. So. Ro) such that
h f(: 0 So.
for some parameter value 0" and the pair (Ko, mn) is assumed known, (In the following, we omit explicit notation the dependence (Ko; mo) for brevity))
Assumption 2. The vectors 0% (8.B) and 1 (8.3 belong to compact subset
(R
C (3,3)
defined by
(3.R) = Pko-1 x B(n. B x Doo x B(p, (81) x B(p, [R)) x Dis-1
where the components are defined as follows:
Pk-1 = [(T TK) E [0, 1|sK E Th = 1] (set of proportions):
B((,r) E R |||||| S 5 with ||>|||
B(p.g.r) {A
|||||| S p}. with norm
||||
sup ||yl|::
VERY, 11yl==1
Dr is the set of positive definite matrices with eigenvalues in [sm, SM] (with 0 < 8m 8
Assumption 3. The optimal parameter 0 Ba) is an interior point of O (S....o)'
The assumptions can be extended the SRUW framework. We define the optimal parameter as 0%. and the maximum likelihood estimator as Ov. In this modified framework, the assumptions become: Assumption 4 (SRUW). There exists a unique (Ko, mo; ro, la, Vo) such that
h = f(; (Ko. To. Noo))
for some 0-, and the model (Ko, mo. To: lo. Vo) is assumed known.
Assumption 5. Vectors 0t and Oy belong to the compact subspace Oly of the compact set Oy, defined by:
0 = Pko-1 x (3) x DKo x (U)) x B(p. R. (U) x D 1j x B(n W)) x D (
Assumption 6. The optimal parameter 0z is an interior point of O'y-
D Proof of Main Results
Throughout the proofs, we denote the full parameter vector of the K K-component Gaussian Mixture Model (GMM) by (л,, {1u ) ) where are the mixing proportions, f4, E RDD are the component means for D variables, and V are the component precision matrices (inverse of covariance matrices Ez). The true parameter vector is denoted by a". The negative single- observation log-likelihool & (yn; a) = - In felust x),, and its sample average is (x(a), The (y:o)
score function components are Sy ax))
do, a=a*. The Fisher Information Matrix (FIM)
I I(a"), and the empirical Hessian is (a). The set of true non-zero penalized parameters in a is So, with cardinality $u- The error vector is A = â a'. For ranking, O (j) is the ranking score for variable j based on grid of regularization parameters Gx. Standard notations like [[] for expectation, ||1,
I1 || I1>o for L1, L2, Lix norms. Other abbreviations include KKT (Karush-Kunn Kuhn-Tucker, RSC (Restricted Strong Convexity), and SRUW (for the variable role
framework).
28

---

D.1 Proof of Theorem
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 G.MM representation for the SRUW model. Second, we apply the identifiability arguments for the complete- data SRUW model similar to [38], which itself relies on the identifiability of the SIR model [37]
Proof. Let M1 with M2 E (K..)))) The complete-den density for an observation Yn under Mj and Equation (12) is given by:
(K, m,r, I, V. 0) and M2 (K", "mmm.,",, 1", V 0*) be two SRUW models
K. K.
Σ ps felus (B:: 2 A ( k,oo):
Analogously, we can derive the same expression for M2. Consider the specific missingness pattern Mo where data are missing for observation n. For this pattern, yo: Yn. 2//00 in. and (n)
Â k.00 Ár. Under the premise of the theorem, if the equality of observed-data likelihoods for
pattern Mo. we get:
K
K'
Σ
Th felust (yn' Dk: As) = Σ Pke felust (yni Ue AR),, Vyn E R. ===1
From theorem in [38], since SRUW are identifiable under complete-data scenario this implies K K", and up to a permutation of component labels:
T = Te.
ix = i
and A = A for all 1,..., ......K. K.
(15)
For any y E R:
Freg y a + y 3,52 1: 71,,, y, )+= fehust ( LUW, y² 52)
3
where S = UUUN W. à = ((a)T. (T)
(with |0]
W x RR), and Ω
0
blockdiag( (Sadiiagg((21)))) I). Here, blockdiag(S1) [2, 1) is the block diagonal matrix created by aligning the input matrices Z2 r. The SRUW density can expressed an equivalent SR model density:
fly; K, m, 0) felust (y²) K, m, xx) felust (30" ;ā ( (
faR-cy (3: K, m, S, R, OsR-oy)
where OgR-ca (a, à _, 52). From Section D.I, the equivalent SR densities for the two models, where the second model is denoted with *) must be equal:
feu--q (y: K, m, S, R. 0ss.ey) fek (y: K", .mmm.. 3n."",K" R"",R""0 0%
By the identifiability of the SIR model [37] and three conditions in Theorem 1, we deduce: K K* m = m,, St, = R", a = a", á = à", / 3r, and 12 52"
Since = S* and S = S supposed that, for contradiction, that U n W + 0. Let U n W. Since j E W, for all q E U the rows of (B).) are 0" Additionally, since j E U* there exists some E E* = R such that (3²) # OT. However, we have established / == 3r. This is contradiction: ())a) must be simultaneously O and non-zero for all E sR. Therefore, U n W 6. By symmetry, Un 8/* 0. Given U W = U U W the disjoint conditions imply C U and U' C U, hence U = U Consequently, W W* This establishes V = V*. The parameter then follow directly from gr. Ω Ω¹, condition in Theorem Thus, the
equalities a = a².B B", 7="", comparing the components of à á" SRUW parameters under MAR are identifiable.
Ω 12',0 r I and = =
Now, extend this to the full model which includes the MNARZ mechanism parameters & {vx & K--1 The full set of parameters for the MNARZ-SRUW model is (0, 사b). Suppose that the two models MA and MB produce the same observed-data likelihood for all observed data (ve, C.):
LMNARE-SRUW(PP Cni @dhta MNARF) (MNAR cn 0 + MNARG) V(yh.c))

---

By Theorem 5, the equality of observed data likelihoods under MNARZ
MANNAR--RUW F Cn: @data A = IMNARA ( you Cn. O
BBB
MNARZ)
for all (yh: C.) implies the equality of the likelihoods for the augmented observed data yn (yh: C.) under the MAR interpretation where f(nn Znnk = 1; 12)) ) is treated as part of the component density:
KA
(Ta)xfa(ull2. =
=1; 1: (x) (baa))
K
fw = 1; (OSR [Cn) = 1; (1b))))
(16)
for all (yh C.). Let ge(y Sai (0A)s, (bA)&)
fa (3. Znn
1; (OSRLW) )x) f
1; (wb)))) This means two finite mixture models on the augmented data (yh, Ca) are identical. If Lk=1 (T ())))=(()==
=1 K2 (T))) (+), and the family of component densities {gx) (and (gke }) identifiable and satisfies certain linear independence conditions ([68), [2]), then K = KB = K, and after permutation of labels, (TA) (TT)) ) and (x)-) = g (-) for all k. Then, we need to en- sure the family of component densities hg (gh, Cni 0Nr 1hp) = ( Pmn = 1; (f.l)(((clll)) = 1; tha) meets such conditions. Hence, we further assume:
1. The parameter sets @SRUW and MNARK are functionally independent
2. The support of yn and C₁ is rich enough such that equalities holding for all (yn: Cn) imply functional equalities. This implies there is sufficient variability of yn and e..
These assumptions hold because v parametrizes only f(c,, (z) and @ only fz hence there is no functional overlap. Applying those assumptions to Equation (16), we conclude that K^ KB K, and up to a permutation of labels:
(TA):: (TB)) for all k,
(17)
and
fa fu Znk 1; (OSRUW) (c)f(ccl2nk 11:
(18)
fe(VElZnk = 11 (ASRUW (O (c))(c =
for all (y: e.) and for each = K. Now we need to deduce that (GSRUW): (OSRUW) ) and (bA)R = (WB)). Equation (18) can be written as:
A Bz = A2(0))
where (1() = fa = 1; (OSRUW))) )s). (c,,) = f(cn = 1; (1b4)). By defined
assumptions, this equality holding for all yn C. implies relationship between these functions. Since y, and cn can vary independently and the parameters (). (v'A)& are functionally independent, J Ar Bi dc AgBady 1, and B₁ A₂ pointwise, and assuming A1, B1, A2, B₂ are non-zero almost everywhere on their support: Pick a cy such that [(c), kk; (vA)) # 0 and f( k;; (bB))) * O. Then for this fixed cf:
fx Piik = 1; (:SOSUU)) )&) CA fx 1; 10 )&) CB
where CA = f(c. kb; (ba))) and CB f(c))k: (1bB)s) Since both fu fa(v.l... are conditional densities, this implies CA CB and thus
fe(vEl 1: (0SRUW) )k) fe 1; (ØSRUW) )&) for all y4
By identifiability of SRUW under MAR, this implies that the structural SRUW parameters (M..TALl VA) must equal (MB, TB, Ip, VB and the parameters (6SRUW), ) (6SRLW jk-
'Since, each component k. the joint component density factorizes as fa (v zz. 1; 0v) f(c. 1ank
1; vh) and both factors are positive on sets of positive measure, pick cp in that common positive-measure measure set; integrating Ai (vf) B) (co) = Az(y.) over yf gives Bi (co) = B2(Co). hence A₁ Az. and then Bi = B2.
30

---

Since CA CB means k: (WA)s) = f(c,, lk; (#B))), and this must hold for all ch (by choosing different fixed c₁ by varying yn first), this implies:
f(en k; (bA))) f(ealk; (bB))) for all ea.
By idenitifiability of MNARz, this implies (vAA) (+B)).
Thus, we have KA KB K, and up a common permutation labels for L..... K:
AT (TT))
(MA,TA, A. VA) = =(mB, TB, LB, VB)
(0XRU)) ) =(0B )r
(VA)K (UB))
This implies the full set of model specifications and parameters for MA and MB is identical. Therefore, the SRUW model under the specified MNARz mechanism is identifiable.
D.2
Proof of Theorem
Proof Sketch.
Suppose that the observations are from a parametric model Po (f(:;): E θ), and assume that the distributions in Pe are dominated by common a-finite measure with respect to which they have probability density functions f(:: 0).
Now, we need to prove the consistency of the sample KL, which will lead to the consistency of the BIC criterion in both the SIR and SRUW models under the MAR mechanism. First, we prove the consistency of the sample KL in the SI model:
Proposition 2. Under Assumption and 2, for all R),
N
h(ya) fVy; 0 ((RR)
P
DKL h, f(; (i.R))
In
Proof. Let O C. (I...., D} denote the indices of the observed components. Recall that the full covariance matrix for the component Ak. is built as:
Eiglls (E
j. S, E S, E S
(ATE))s
E S I S
( A E
j, S
For indices j. I S: As before,
Ak.it = Ex.it-
Thus, eigenvalues of A are in [sm, sAs]-
For E S, E 8c:
A = (ExA)st
The norm of the product satisfies
Σg S EAll -
Since Σσ Dj8)+
Σx Xmax (E+) SM.
Given B E B(p, (R,
A S.p. p.
Thus,
A S SMP.
This bounds the norm of the (S 89) block of A
31

---

For j.7 E 8%
A = (S + A Dzt
We can bound the eigenvalues of this sum using the spectral norm:
|| + A A A |2 + E
Since Z E Dist)
[2 Amax (S2) S SM.
Next,
|AT S |A"
| ||||| | ||E S sas...
Thus,
Ω A E S SM + pMSM = AM (1 + p²))
Given the above bounds, there exist refined constants s. > and
M = m p²)): SMP. sM} = SM (1 + p²)
such that
m k SML.
Note that AM will depend on the sizes of the blocks, and it finite being dependent on 8Ms P. and the structure of Ω. B, and Σg.
Since 0, is principal submatrix of A and A symmetric matrix, by the interlacing theorem:
8m Amia k,o0) Amax k,oo) I $31-
Thus,
8m T|a|
A k:00 I 8as
and
k
/101-
8m
We are now moving to the main proof. By hypothesis, the true observed-data density is given by
h(V")
Σ
The 1 y Mkoor Ah.00)
1=1
and if we show that
Ey [1 00.
then by the LLN,
Σ In((v))) P By [In h(y))
Moreover, the observed-data likelihood under the model is
0v Σ mh & ( Ak,00)
and we later show in Proposition 3 that
N
In(f(s.:0
P Eye [In fay (( (G.R))] n-c
To do so, we must verify that the class of functions
F(S..) = {1r + In (): O (S.R))
32

---

satisfies the conditions for uniform convergence under Assumption 2. In particular, by Assumption 2, the parameter space O (8,3) is compact, and for every y E R [0] the mapping
+ In ()
is continuous. Next, we verify that there is an h-integrable envelope function F E F(S.R). Recalling that for y E R (0,, the Gaussian component density is
1 VE:o Ak,00) = (2xr
Vk,0:
[0]
A R
exp (- (1) -VK,O) k,00 10000 (v" vk:0))
We can derive the bounds for this as follows. For the upper bound, since Amin k,,o) 8m
101 4 A kool 1 S (3m)
With exp(-.) S and A,,0 is positive definite
+ y Uk,0, Ak,00)
10 S (2x)
(8m)
(8mm)
191
Summing over l and using Σk=1 Th 1,
101
s(y::) S (2Tmm
Taking logarithms:
0
I In fly";
In (
2
For the lower bound, for each k;
In (n Vk,o. An.00)
10
In
1
In Akoool
1
(v00 -Vk,o)" "' A
k,00 (3. - -k,o).
2
2
Since Amax( A kk,0o) I SM.
In A kooo S (n)m)(
and A 1000 1 Fm
11
1 Vk.ol1
vkoo) s (y) 1k,0)² " Ar
5m
Using |/ - ukkoll? < 2111)/112 + 1113)
2(118/112 12+ + || 112)
-k. T² A 1.00 (3
vkoo) -
5m
Given the construction of the mean vector uy for each component A:
pkg
j
Vkj
(a+ ( Ax))
if E S";
and knowing that belongs to the closed ball B(p, 1, D set of |S matrices with norm bounded by p-we seek to derive a uniform bound for uı. For indices j E S, 1k₁ Pkj, and since the parameter space Or is compact, the cluster means μky are bounded by some constant 17 > For indices j E S we have Vıj = (a + PxA), = a; + (mn. Using the elementary inequality (u + v)2 < 2u2² + v²), it follows that
|| < (10m]1+
Given that a lies in B(p, 1, (S)), ||,,11 < p² for each j E S addition, the term AHaA) can be bounded by [Mell2l1 A 12, where A, is the j-th column of A. With μ4 belonging compact set B(n, [S]), || 1 n7 and because A is derived from bounded parameters (including
33

---

3 with norm being bounded by p²), each column A, is bounded in norm by p². Consequently, (1xA)yl & p² + +ppp² p²(1 (1+ + ) for each j E S Combining these results, for any j E ES,
|| S (1"))))
Thus, all entries of vs, irrespective of whether they correspond to indices in S or S", are bounded by a constant that depends on y, p. Hence,
|1 Σ ||xk
3 Σ 2p +1²)
(1 +1/)
uniformly for all k. Consider again the Gaussian density for the observed variables
In /b y": y Vk Ag,0n)
0
ln(2
1 In
A R,oo|
)k A
k (2 - Vk.O)=
2
2
2
To control the quadratic term in the exponent, we note that
11 Vn.oll (
Given that | = EjEollul Uky 12 2|00pp +22 is uniformly bounded, follows that
11 vn.01 < 1100
This bound ensures that the quadratic form (y- -k,o uniform boundedness of vı thus allows to assert the existence of integrable envelope function F(y²) that uniformly bounds In ((y 0)] for all 0 within the compact parameter space Oy, Thus.
K,00 (y) Vk.o) remains finite. The
1 21 21 n)))
In ① ( V Ak,00) 2
In (2()
In(SM)
2
2
2
8m
Simplifying.
10
|| 2|0) 2/0)p2 17))
In 2 y:: uk.or A Using Jensen's 's inequality over the mixture:
1n(275M))
2
8m
[0|
11 22/00p2 +
In f
In(2
In
8m
Since Σk Th = 1.
|0]
1 22/+220pppp00 + 1?))
In f(y":0) 2 In 1n(2T3m)
3m
Combining the refined upper and lower bounds:
In(2TT8r)= || + 2/0 p² + n?)) S (('";
(0 In (2
8m
As final note, these bounds also rely on the eigenvalue constraints on Σ,, Ω. and the norm constraint on B.
34

---

Now we prove that the envelope function F, which is related to the upper and lower bounds above, is fi-integrable. The true density hzy) corresponds to a Gaussian mixture model fly: (s),Ht)) with parameters in compact set. The observed-data density for y", obtained by marginalizing over the missing components, is given by
K
h(y") = Σ Th 1 y: (
Sk.00)
kml
To verify that the envelope function h-inttegrable we need show
/ |1"1)>>(">) dy < co.
We proceed by examining the second moment of the observed components. First, observe that
/ 11 dy Σ Th
/ |1801122 1b
( Vk 2k,00 dy".
k=1
S 2 2 II2 + + 2tr
by using Lemma. By the compactness of the parameter space Og and since Ek= Th and the bounds derived earlier, we have
B ( dy 101pp( + 28M S0, 00
Therefore, F is h-integrable, i.e..
J ((g)))(()))) dy 0o,
since f (1'") dy 00. Finally, because [n(h()))) F(y"), implies
Bl|in(h(y))] / (n((())))))))))))))))))))))))))) dy s F(y")h(y")) dy <0.
Hence, the envelope function F is h-integrable and E|| In(h < oc. Thus, we can apply the law of large numbers to obtain the consistency of the sample KL divergence.
Proposition 3. Assume that
(yf...., yh) are i.i.d. observed vectors with unknown density h.
2. ④ is compact metric space.
3. 0 E θ > In(j 0)) is continuous for every y E IR a1,
4. F is envelope function of F (In[f(:))]: 0 @) which is h-integrable.
5. 0* argmaX@ge DKL hh, f(: 0)]
6. θ = argmax Σ In f(y.: 0).
Then, as n 20,
N Σ In f(y%: (s,)) E Ey In f(y": (3.R))
Proof. We consider the following inequality:
N
E In [(y"
(G.R))
In fly: (s,)
35

---

< In fly,; (s.R))
E LI f(y";
f 0 ((3)))))))))) )mm)))
+ sup
In fyy - Σ
N
In f 0)
aa.
By the definition of 0 (8.R)* have
E ln fly: (i.RY)
13 In f(y": 0 ((.x)) 20.
Note that,
N
E In fzy (,,,)) = E In fly"; 0 (sn)) E In f3y ((.R))]
In f 0 (3.R)
N
n=1
N
In fly (sR))
Σ m=1
In f (8,R)
In I 0 (BR)) E In (3y 0g)))))))))))) )
For term 11, sine maximizes the empirical log-likelihood,
N 0(..()
1 N
Σ In(S(Xa (o,R))
which implies that // S Therefore,
E E[In f(y;; (ARR)) E [n f(y: : [smm
s + ///7/-
Moreover, since both | and |111| are bounded by the uniform deviation:
| < sup
In fz (3::0) 0)
In f(v: 0)
8 (5,R)
N n=1
and
17 < sup
In f 0)
N Σ
In fly 0)
0c (
Hence,
E In fly
0( ((RR))
- -E In f 0 ((.R1))
$ [1| IIII
1 N
< 2 sup
In fn 0)]
Σ
In f(y 0)
(5.8)
Hence, the left-hand side bounded by three times the uniform deviation:
N
N
B 1as 0 (B.R)
1
Σ In fly::
0 (s,R))
< sup
E In f 0)
1 N
Σ
In fw. 0)
N
n=1
By the argument in the Proposition 2-namely, using the compactness of θ (8,R). the continuity of 0 + In {:; 0), and the existence of an h-integrable envelope F-the class F(B.3) is P-Glivenko- Cantelli by applying Example 19.8 van der Vaart (1998) to conclude on the finiteness of bracketing numbers of I under the assumptions. In particular,
N
sup
In 1()
C 0
(
(8.3)
==1
22000
36

---

Therefore,
N
N In f(un 0 (s.R))
Σ
& n-
E fn
0 (s.)))
==1
which concludes the proof Proposition
Now, prove the consistency of the sample KL in the SRUW model: Proposition 4. Under A Assumption 4 and 5, for all V,
1
In
h(gg.)) f(y:: 0y)
P.
DRL h ((:
n=1
Proof. We carry the proof in the same manner in the SR model. Let V = (S. RR, U, W) and O (1...., ;) D} denotes the observed variables. We still want to apply the Proposition 3 with the new family
F(m) = (Im]((10)) o E 0
similarly to the proof of Proposition achieve
N
In(f()):
& In 0.)]. n-xc
N
To do so, we must initially verify that the class of functions
F(v) y" In f(y/ 0
satisfies the conditions for uniform convergence under Assumption 5. In particular, by Assumption the parameter space O₂ is compact, and for every y E R|0l the mapping
f - In ()
is continuous. Next, we verify that there is an h-integrable envelope function F E F(v). Recalling that for yd E R [0) and given component k. the current Gaussian component density is
T ( As,00)
(2xr
1011
exp
(v iik,,)) À k,oo (y) ik.o))
We will bound the density function usual. From Proposition 2, we know that there exists 8m > 0 and 3y = max(sM(1 + (²), SM) = 8 (1 + p²) such that
5m A k SMLL
Thus,
Sm Idll 1
10 I ā I|li
and
A (A,0o)
8m T|o1-
Given that E D/w and is independent of any relevant variables, we have that the principal sub-matrix I is bounded by sa
Amax (Too) S SM-
Since Á 8,00 block-diagonal, its eigenvalues are the union of the eigenvalues of A koo and Too. Hence,
Amax A k.o max(sMM( + p²), SMn (1 p3).
We define the upper bound by
8M = (1 + p²)-
The lower bound of the structured block, together with the lower bound on the independent block, implies that the covariance Å k. satisfies
8m I < Å K000 3 3 I.
37

---

Because -1
I 1. I, it follows that 2m
k,00
1180 - ol
(y - 1kkoo)
k,00 (3 -ik,o) S
5m
Using the elementary inequality
11 S
110x:ol1)) +1
we obtain
2 2(118112 1k.
(y² vk,0)* -kko)
2 A
K,00 (30 ixk00) 2
8m
Now we process the bound for Vk.o- Recalling that if we denote Day denote the number of coordinates in the structured part (i.e. U U); then
|k
11
Σ |kall s p
JED
uniformly bounded and l/al1 & 172 since y belongs to the compact set B()), W[), it follows that
||0x,ol1 s Dsu,op + +n) +
We deduce that
Dau, 112) 1/2
In ① ( ix..0 k,,00)
In In(2T 5m))
5m
Using Jensen's 's inequality over the mixture and
Σ k= Th 1, we have
K
10
1 Dau o2 (1 + 172) + 172
Tk
In 1n(2TTM)
2
3m
k=
0
| + Dis.op + 172) +12
In 1n(2TTm))
2
8m
Therefore, each family's member in F(v) is bounded by
0
0 Dsu Dsu,op
+
|0 In ln(2rmm)
In(2T5m)
In f
8m
Now we prove that F is h-integrable. The true observed-data density under the SRUW model is given by
h(y)) Σ me @ ( ok.
where the "*" indicates that the parameters are the true ones and the density is that of a Gaussian mixture with parameters in compact set. In our derivation we have shown that, for any O in the family F(y). the log-density satisfies
[ (2(
|| Dov p² + 1/2) + 112
10 S f(y; ; 0) 2
In 6m)
Sm
In other words, every function f(y"; 0) bounded absolute value by
F(y") In (2 + ||2 + DEJ.o p² + 112) + 1/²
2
8m
Hence, for all 0 in the compact parameter space,
In s S F())).
38

---

To verify that F is h-integrable, we must show that
f F(y") h(y) dy <0.
Since the envelope function F is of the form
1
Fly") Co
8m
with
[0 Co In In(2TTM)) + Dsu.o p² (1+ + 1/2
8m
it is clear that
/ 1 h(y)) dy 8m
Fly") dy S
Because hyy") is Gaussian mixture with parameters in compact set, standard properties of Gaussian mixtures guarantee that the second moment is finite, that is,
f |1s 1(y)))) dy co.
Thus,
/ F(y") dy <0. Finally, since for every y we have
[nn(yg)) s FF())
follows that
E [[nh(g))) = / h())) dyy & f F(y))((y)) dy <c.
Therefore, the envelope function F' is h-integrable and E[] In(h < co, and consequently, we 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 = (8, R, W) by
BIC(V) 22 (0v) Ecyy log N,
where Evv) the number of free parameters. Let Vo be the true partition and set ABIC(V) BIC(Vo) - BIC(V).
Let V1 V # Vo Dku(h, f(: f(; 0}. The Identifiability Theorem 1 implies every V Vo belongs to V1 For V1,
[(yh: [(oRiBb))) h(ya))
)
ABIC(V) 2N
Σ In
In
[Eer) -(vo)] log N
N
h(y))
N ==1
n=1
(19)
To prove the theorem, we will prove that:
VV e P/ABIIC() < 0] + 0
N-100
[[yg h(y2)
Denoting (V) DN In
.M(V) - DKL, [(; 07)), from Equation (19),
we get:
P(ABIC(V) <0)
= P(2N (Vo) - MM. + [Ev) - [(v)) log N <))
-P(MN(VV) M(V.) + .(V.) M(V) + M(V) Mv B(y) - E(vo) log N
< <0)
2N
[[v) -(v.)]
log
<P(MN((V) -M(V)) > c) PPM(V) M. > - -MV) +
2c)
2N
39

---

From Proposition 4, we know that
1
N
In
M(yh)
Pg n-000
DEL h, f(: ((:
N
f(y;; 0v)
This leads to My prove the theorem.
P n=000
M(V. Similar to the proof for the SRUW model in Maugis [38], we
D.3 Proof of Theorem 3
Throughout this proof, we analyze the estimator defined by minimizing the negative log oggliikeihood penalized. Accordingly, although we denote by
a(y: in ay)
the usual log-likelihood of single observation, consistently work with its negative -f as the optimization objective. To streamline notation, we therefore use the term "score" to mean the gradient of the negative oog-likelihood:
(y a) de(ynia) da,
This convention differs by a minus sign from the standard statistical score, but it ensures that all gradients, Hessians, and Fisher information matrices below taken with respect the minimized objective -Cn-
We start with presenting establish fundamental properties of the penalized GMM estimators and the variable ranking procedure.
useful lemmas before proving main theorem. These lemmas
m- Tm T (yn's Mm: Em, be the p.d.f of K -ccomponent GMM for an observation yos where (T,, (um)k=1 {Ex ( )n) represents the GMM parameters. Let i Σx be the precision matrix component The soglikelihhod for observation yn is (yyn a) = In fetust (aasr(()r (yni c)). The responsibility for component k given observation yn and parameters or is la Lg(yni (r)
Lemma (Score Function Components). Let (oa ax) Σ K
Te Teynsti Es), The score components, Sy(yn; O(")))")) where oz Jawel (nii))
a parameter in and a. are the true parameter values, are given as follows:
1. Mixing proportions Thž When parameterizing म = EØ= for Tyy [l,...,,K K 1}, is:
K-1 Ti.
the score component
(e(yn;("))
Sr. (Yn:
Th
2. Mean parameters (px)a: Let (Hua)d be the d-th component of the mean vector pę- The score component is:
S(pa) a²) (yni a")(" - =a))
3. Precision matrix parameters (wx) Let (Wa)re be the element (r, 8) of the symmetric precision matrix R. The score component is:
(2 drs)
S(ve).. (Yo; ai") t/ (y); a²)
((E)rs (yn pp)r(y - p%)s)
2
where dra is the Kronecker delta.
Proof of Lemma 1. Let (Wni yn is Elynn ar) In Emm Tm 1 (yni ox) Syyn::" da; a=a".
(N a D(yn' Mm: Em). The log- likelihood for single observation
The score component for generic parameter ay is
06(():a))
40

---

Score components for mixing proportions k E (I....,K 1}:
We parameterize T = - ΣØ=1 Ty. For
allyn:a)
K
K-1
Σ
Tm Pm (Yn + 1- Σ Tj
--11111 j=l
(x(Unia)
ŌTk
felsa (yni a) OTK
m=1
D (1) ax)
T (Yni
Seluse (Yni (c)
Since tm (Y:: a)
Tue Pm (mm (W..ic)) we have
om (yria)) TamlVanaa
fm (M) Thus,
Jomn (B)
)
4 (yni
(xyynj ox)
ŌTk
Tkk
TK
The score component at o² is:
4a va(unoa)) =0=(""("00)0)00)) 1x (yn3)())))
L. (ynn:)))
(x()n:a"))
ST. (Yni (=)===)==>
oi
K
T
Score components for mean parameters (ux)a. Let (pu)) be the d-th component pk.
a
K
Σ
Tm Pm (ymi or)
0(z))
Scluaa (yni a))
a(pss)a
Tk
(
Sclual (Yni (z)
a(ps:)a
Since 00x 0 ts. and for 42. In (3 a.) -C -2 -M0)T T p4x), have:
4b a(ux)a
8(sx)a
b In t2 ()
(V mi))la
O(uxa
Substituting this back:
df(y); TA T T x
(V px))a 4x (g): =x))a
H(u))a
fetual ()):
The score component at ax is:
S(nw) (yn; -tx - 14))da
Score components for precision matrix parameters (1rr)s Let (Vrs be an element of the symmetric precision matrix L
de(yn:a) x( Selust (yn
PDe (1 ar) 1 k
1 Oln (a((naa
T a(vr)).s
T
The log-density of single Gaussian component is In 2x ar) = Cy + 1 del ( p&) 1 p.).. For a symmetric matrix R the derivative with respect to an element (Lx)ss (using the "symmetric" derivative convention where ax... r * s) is:
means varying X,s and X,r simultaneously
i In det VT)
(2 (v "rs =(2 (22
V Ire
oly px))
( Ma) (2 6r.)(yn (k)r px)
0 1
Therefore,
i In (2()n:a)
= (2 (0,--,) -8,)) 7 6,.) (x)r(2n - px)s
0 k)r
41

---

So,
ol(yn: to (3n x)
(2 (2-5rs) 8rs)
((Ekrr - (Yn - (3)) - (4)a))
0 krss
2
The score component at ox is:
(2-ôra)
(yn Ix (y(nnn"))
((Di)rs - (yn - -1b))(( - p()))))
2
For the two lemmas below, we will use the following assumptions:
Assumption 7 (Identifiability Smoothness). The GMM density (elust oi) is identifiable. (y: oO) is three times continuously differentiable w.r.t. in open ball [B(x'' ro) around o².
Assumption 8 (Compact Parameter Space & True Parameter Properties). af is interior point of compact set O' C B((a"rro). This implies:
Mixing proportions: z 2 Tmin > for all 1...., for some constant Tmin (0, k.
Means: |eella S do for all k.
Covariance and Precision Matrices: The covariance matrices ΣS have eigenvalues X(E)) such that < 8m I X(E)) S 8M < cc. Consequently, for the precision matrices
V = (>i)) their eigenvalues (VV))) satisfy 0 < 1/SM S > A( I 1/8m < co.
define min Sm. pmax SM. And for precision matrices, 9mnn 1/sM, mmax = 1/sm-
All Ol satisfy these bounds. Let Ly be upper bound on the norm of the third derivative tensor of (yy a for a Oz such that [y)) co.
Assumption 9 (Data Distribution). Observations y are i.i.d. from felust (y:a²)) O²), For each yo. there exists a latent class variable 7, (1,..., K} with P(z, k) The such that Yn 1. = k ~ (N(, D
Assumption 10 (Bounded Posteriors). For a E Oy the posterior probabilities ( = T.o(y/me... E.) satisfy 0 < (x) < E (,0(s)),
The empirical Hessian is (a)) = (a)) The Fisher Information Matrix is I(a") Ea (Va-(V o.). Moreover, we have the following definition of the true support of penalized
parameters:
Definition 3 (True Support of Penalized Parameters So). Let a = (TT, [ug (1i)k be the true GMM parameter vector. Consider the penalty
K
K
P(a) = = P
1
+
k=1
k=1
We define the true support set So as the collection of indices corresponding nonzero parameters in a" that are subject to penalization:
Sp Sy
US
80
Sol
s
Here:
1. Support of means
s = 15 15d5p, (1t)2 + 20, 82 -
2. Support of precision matrices
S (kkrr,s) 15 15 (Wi)rs 40.
si = =15
42

---

Assumption 11 (Restricted Eigenvalue Condition). For any parameter increment vector A indexed compatibly with a, we Ass So for its restriction to indices Sg and A (So)e for its complement. For 2 1, define the cone
C(co. 5o)= { |1 (56) I coll A ]
Assume the FIM I(x²) is positive definite, then there exists K1 > 0 such that, for all A E C(co,
A A(())A KM A
Remark 1. The constant is chosen to match the cone where the estimation error â - a lies. A sufficient choice co Z As+1 when the regularization level satisfies X Z Ao VIVex (a"))l with Ao 1.
Ag
Assumption 12 (Uniform Hessian Concentration). et da (K - KD KD(D + 1)/2 denote the number of free parameters in a. There exist constants Cn > 0, 6 Jf such that, for N 2 80 In da- with probability at least ct d
> 0 and radius on > 0
8o In da
sup aEB(a' r AEC(co.So), ||
A (N - A S CH
The constant Cи depends on the bounds in Assumplion and on moment bounds for derivatives of E.
The radius can be taken as öR x / In da/N..
Lemma 2 (Gradient Bound). Let a = (w (i. pk Vi..... xx) be the true parameter vector of K componentt GMM with D-dimensional components. Suppose Assumptions 7-10 hold. Let da be defined Assumption 12. Then there exist absolute constants C₁, C2 > 0 and Cg such that, for any N 2 1,
Inda N
P [VVen(a')l S SCaV
2 C1
where (x(x)
a) is the empirical negative log-likelihood (recall lly; a) In felus (y: a.)). An explicit choice is Cg- /2((2 + 1)/co Vmax where co an absolute constant in Bernstein's inequality (e.g., cu = 1/2). This bound holds provided 2 CB ln da. with
2(2 + 1) 2
CB co ming:: by 70 ( (B(S)) 21
2 maax max, E[ST,) and by is the sub sub-exponential scale of Sjn (by = 0 for sub- Gaussinn compo- nents).
Proof of Lemma 2. By Assumptions 7 and 9 (regularity and correct specification),
Ex- (s) a²)]
for all j L..... da
Let v3 = (SS)((n: a ) be the variance of the score evaluated a... We bound variances and tail parameters for each parameter block.
Mixing proportions Sm,, (yn; a²). By Lemma 1,
Sz, (yni (²" ax)) txx "yn::: a²)
tx (1):: a*)
By Assumption 10, 0 4m (ynyion a²) 1, and by Assumption 8, nm > Tmin
0; hence
SSp (Bn a<)] < 2/ Tmin. and byg 0.
Thus Sny is bounded sub- Sab-Gausiinn with 2 = E mp
< (2/Tmin)
Mean parameters Sim (Yni a*). From Lemma 1,
S(w))a a") - -ta a²) (Wi mil)
43

---

Let Unkd == (Vi -1x)) Conditionally on 2n = l, we have
p4 tEe with ~ No,, E;}.
Then
Var 2n = = ed D Rea I | Eilla
O muax.
and
UnUdd In =n== led ( - p - 1 l - Hilll I Omnx (2n)-
Hence BlU nkal S 6 + 4n Since |t -
Vin B [ (m
< Rual S
max maxx + 4n7²).
Moreover, Unkd is (conditionally) Gaussian and hence sub- Sabs-aass;;n; with the bounded multiplier ts. S(pn)a is sub-Gaussian, so b(r))a 0.
Precision entries (( a"). From Lemma 1,
S(ww).. (yn: a²) - (1xxix)))) (2 (2-8re) (Ei)re (yn (1) mi))))
Let Vnnks= = (Egrrs (yn (1))(un - M%)s- Each (yn 14)), is sub-Gaussian with parameter S max +77 so the product (yn (4)r :1))s s is sub-exponentill hence Vukrs is sub-exponential.
Since [t2) 5 1, S(wrar is sub-exponential. Thus there exist constants Cv. Cy.b >> 0 such that
((w.). E
S [(wwr.) S Cv.y (o n/2)2
Cv.b (o max n7²)-
Union bound. Let viax max, E[S and buaxx = max, {b, : by * 0}: by the above bounds these depend only on (Tmin: 77, 8m: 835). ). By Bernstein's inequality, for each coordinate j and any >
N
Sjn 26 < 2exp (-coN min
P
by
n=1
with the convention that if by = 0 then min(-) = f/v33. Cg In da so that tn I minj;b, #0 v2 /bg. Then for all j.
Set ty
Cg C
Ind. N
and assume N 2
N Σ
Cu.Ntyy 2v²
Nt
Co 2
P
Syn
2 tv
I 2 exxp
S 2exp
de
2n
n=1
union bound over j ,,..., de yields
P (( 2 (N)
2 2da
2n
Choosing C₁ and C, so that ca = C2 - gives the claim, i.e., C, = 22((2 1)/c0 Vmax 2
Finally, the condition N 2 C3 In dg is guaranteed by taking
C2
2(C₂ + 1) D
CB (miny:
tniny:,,z40 2, 13/b,)²
co min b, 2o (Eb)),
Lemma (Parameter Consistency for Penalized GMM Estimator). Let & be any local minimizer of the penalized negative average log-likelihood
Q(a) = (x(a) + P(a),
where
N
1
(a(
In fetuse (Un: ox),
P(a)
Σ
p& 11 + +p
Velli.
N
k=1
44

---

Let a' be the true GMM parameter vector, and So be the true support of the penalized parameters in a", with sparsity So |Sp)- Let da be the total number of free parameters in a as defined in Assumptions 12.
Suppose Assumptions 7-12 hold and the gradient bound in Lemma 2 is satisfied. Choose the regularization parameters X = p Xchosem where
Inda N
Achosen 4o
for constant Ao > (e.g. Ao 3). Then, provided N 2 In day with probability at least R
1 Cqrud da 2 c da ?
1. (Lg-norm consistency):
2(Ao + 1)Cg
80 dn
|a - a0112 1a-a0l112
K1
N
2. (I--norm consistency):
4Ao(A0 1)Cg
In da
la ---al a=1l1 S - 1))s s0
N
Here Cg is the gradient bound constant from Lemma 2. and 14 is the restricted eigenvalue constant from Assumption 11. The factor Ag is user-chosen constant controlling the regularization strength.
Proof of Lemma 3. Let ā be minimizer of the penalized negative log- log-likelhhood Q(a:) = (x) P(a), where P(a) Ex | | and (x( N/ E E(yni a:) with E(y: a))
- In felua (y: a)). Since o minimizes Q, we have
(x(@)) - (x(x') S <P(a") -P(a).
(20)
Let A = - ax and assume â E Oy so that a A E Oy. A second-order Taylor expansion gives
(x) - Cx(aa) (VEN(a"),
A HN + L 2
for some to E (0, 1). Write a² + A. We lower bound
A ATH à A by decomposing
2
A Hx A
A (T(@) I(ax') A H (ā) I(ā)) A.
From this point, all bounds are stated for arbitrary the cone C(co, Sn); the fact that the actual enror ā - o² lies in C(oo, So) will be established later by Lemma
By Assumption 11, for A E C(Co,So),
A A(a" A 2 K3 11
By Assumption 7 and compactness (Assumption 8), a + I(a) = E(V? ) is Lipschitz on B(o", ro) n Og in spectral norm: there exists L, such that
[T(@)) I(a 1)
Lılāa a" Litol|Alk2
S Lu A
and therefore
A (I()) - I(aa))))) < [I(@) I(a 1|2 A < |
L A E
provided
A S ÒR. By Assumption 12, if 22 S ôg and A C(co, So), then with probability
least
1 cf d,
In d
A (HN(()-I(a))) A Cuy
A
N
45

---

Combining the three representations, for A E C(co, So) with | S 8R.
1 2 R1 11||a|1 Cuy
In da
A HN
||
N
Choose ōn and N so that Lndn S Kr/4 and Ch V Then, with the same probability,
so Ind
N S 141/4(e8, N (n
A (N))
R1 A
Hence obtain the restricted strong convexity (RSC) inequality the cone:
(x(x' + A) - (x((x) - (Vex(a)),
K1 for all E C(co, So), || S 6g.
2
(21)
where K7, - 2 0.
CF da C
Stochastic term. By Lemma 2, if N 2 Cn do, then with probability at least
In
[VEx(a") s Csp
N
By Hölder's inequality (<r. <= 11lloolllll1)
In In N
(Vln(a") A) < [Vex(a |All S C.
A|l.
(22)
Penalty difference. Let A An = ft - 2 and v. Vx Vi. and denote by S and S the true supports (Definition 3). Using the L Li-triangle inequality on supports: for any vector a, b and support
S of a: ||a||] [|b]]] I |las ball1 [|a bay 1I + (0) Since asz) 0:
|| 11 < || A Sp. 11 pm.3splll1
| sell2 I
| 111
Summing over k and writing Sg = 8% US
P(a*) - P(@) S A(1 ll
| (s))-=11))) + p(1 Av.l. 11 (8,)+.))
(23)
Combining, With probability at least 1 - - de C c do combining Equation (21), Equa- tion (22), and Equation (23) with Equation (20) yields, for all A C(c,, with | ||2 S ōn.
KL 2
< C.
In da N
All + A(I
|A (s.))=1)))) + p(1 s.lli - || ,r (8,) III)
(24)
Invoking the inequality Equation (29) from Lemma (the cone lemma, proved via KKT and the gradient bound, thus independent of RSC),
KL 2
A
<
+ Ao
Σ (1
1
Σ XرMم,1-
Au
JES
jE(S)
Since the second term is non-positive (as Ao > 1 and AyD 2 0), we can drop it to get upper bound:
KI.
1
Σ >
Ao
JE80
Taking the common regularization level Achusen = p Ao
Ind with Ao > and note that N
Ejes,
Achosen As, /12 obtain
KL 2
|
(1
Xethnsen Asoll2-
Ao
46

---

Using the Cauchy-Schwarc inequality
A Aso S 50
As l1,, and since
/As,l2
S
||
1|:
LL
22 Achosco
1 +
v 90
2
2
Ao
||||] 2 * we can divide by
A 12
2Actusen (1
1
1
A
12
V50.
KL
Ao
Substituting Achosen A
N
a a²12 5 2A0C' (1 11
1
In da
/s0
KL
Au
2(A0 + 1)C²
so In da
KL
This proves the Lg-rate with constant CL2
2(Ao+1)CC
For the ratte Lemma 4 yields the cone bound A (8,)=11 S Ceone A I1 with Cooe (decomposability via [44]). Hence, assuming all X3 for penalized components are Achosen):
Ao+1 Au-1
A 1 A Sp 11 + A (S,) ll
< A So |l + Ceone
Ao || = (1 + Camee) A1 l|
< 11 Conee) / 90
lla
(by Cauchy- )chhhy))chharz)
< Coone) v 80
2-
Substituting Coone
Aot1. Ao-1
Ao +1 An -
Ao +
Ao +
2.A
+ Coone
1
Ao -1
Ap
So.
2AD
||||
V8ol A 1/22
Ao -
Now substitute the bound for |
2Ao
2(AA + 1)C,
In da N
|a - a''11
50
Ao
KL
4Ao(Ao + 1)C
In d N
sol
(Au - 1)NL
In da N
-CLi(Ao: Cg. KL)
where Czi C.. kz) and La norms.
41o(Ao+1CC, This concludes the proof of parameter consistency in L (Ao-1)Ns.
Remark 2.
The radius ōn governing the RSC inequality Equation (21) is critical and is typically of
the same order as the target statistical error. If || 1|2 exceeds &n. the Lipschitz remainder L A 1 can dominate, and/or the uniform Hessian concentration in Assumption 12 may fail such a large neighborhood of o²" In general M-estimation analyses (e.g., [44]), one often proves an RSC with tolerance term,
(x(a" -x))
KL 1
In de
A
TL
N
valid on a larger set. Under Assumption 12, we obtain sufficienlly strong control to dispense with this tolerance and derive the cleaner quadratic curvature bound Equation (21). Assumption 12
weakened (e.g., only yielding bound of the form CHy/ /n |
Vso |||| A AHH I)
the cone), then a tolerance term proportional to A i would naturally appear in the RSC.
47

---

Lemma 4 (Cone Condition for L1 L--Penalized M- M-Estimators), Let & be any local minimizer of
Q(a) (x) + P(a), where (x) is differentiable loss function and the penalty P(a) is sum of component-wise In penalties:
da
P(a) =
==1
Let So be the true support of a'. Let A ā a.. Suppose the regularization parameters chosen such that for some constant Ag > 1:
> Aol Vol[Ve((a")) Mbl for all (So)-
(25)
and Xj of similar order for j E So. For simplicity, we often set X; = Achosen for all penalized components, where Achosen Z Aul VLN (ax)) If the RSC condition from Equation (21) holds for
Ex(a* + A) Ex(a") (VEN(NV"
KL A
2
then, for Ao > the error vector A satisfies the cone condition:
MH4ر S
Au +1
Σ X/4ر
(26)
Ao
JE(B.)
JE8o
all Xi for penalized components are equal Achosen, this simplifies
Ao 1
(Bo)= Ipen,
A Ippm,,1111
Ao
where || -llpnn1 refers the Li norm over the components that are actually penalized. If all components were penalized with the same xg. this would be (s,)=111 S Act1
Asulli.
Au
Remark 3. In our case, X; > for mean components jldd, and X; = for off-diagonal precision components (Wx)rs and X 0 for parameters not penalized like proportions
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 support. The key is that the regularization parameter for the "noise" variables roff-support) must be sufficiently larger than the corresponding component of the score vector, allowing the penalty effectively shrink noise components.
Proof of Lemma 4. Let P C (1....,,d. d denote the index set of penalized coordinates and write ||az||pen11 Djer (zzl- In what follows, sums and L -norms taken over P. Assume the tuning condition holds all penalized coordinates:
Xj 2 Ag [VEn(a")]]
VjeP, An
> 1,
(27)
equivalently Aehosen 2 V (x") Mloo.pen when a common level is used. 2 Since â is local minimizer of (a)). ). satisfies the basic optimality inequality
(ā) - (x(a") - P(a") - P(ā). Invoking the restricted curvature bound (RSC) Equation (21) for A â α²
x (ā) - Ex 2 (Vex(a"), A) + KI.
2
2
Note that, only the nonnegativity of the quadratic remainder is needed for the cone; the explicit KL > 0 is used later for Latres Therefore.
(Vlx(a"), A)
AA S P(a*) - P(@).
(28)
For the penalty difference, write P(a) = ΣSEp dylaşl and let Sp S P denote the true support of the penalized coordinates. Using the standard L, support inequality with a² and ā = a" + A,
||||| - |bbl = |as |l bs, || | |/1 < ||as, bs. - | ll..
II some coordinates unpenalized, they excluded from and do not enter the cone inequality.
48

---

and since (a") = 0 on P. we obtain
P(a')-P(aj - S Σ
E ト
JE3p
JE(Se)
Substitute this into Equation (28) and bound the linear term by the triangle inequality over P:
(VEv(a"). A) S
Σ
[[VEn(a"))1AA)
&
Σ
14
Ao
JEP
JEP
where the last step uses Equation (27). We obtain
Σ JESo
A
A
Σ jE(Sa)
)
برم
Au
Ao
Equivalently,
(1
Σ
ا (1-
E
(29)
An
Ao
jEBn
:E(B0)
Since the left-hand side is nonnegative, follows that
(1- Ao
Σ
الاي $ (1+ to
E
)E(S,)
jE80
and for Ao > this yields the cone inequality
EMAI
<
Σ x Ao JE8o
Ao +
JE(S,)
all penalized coordinates share common level Achosen- this becomes
|| (Bo) pen.1
Ag+1
As Ippe,,11,, and if every coordinate is penalized equally it reduces to A A
Ag-
As,l1..
Ag
With the cone condition Lemma 4 established on the penalized coordinates, we now consolidate the standing assumptions for 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.
analysis (Lemma 5) and for Theorem 3. The
Remark (Alignment with SRUW assumptions). Assumptions 4-6 were introduced for the SRUW model in the MNARz setting. 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 (Ko: mo. To, lo, Vo) explicitly here because no model selection over SRUW structures is performed in Theorem
Compactness and interiority (SRUW-5-6). These correspond directly to Assumption 8 below (compact parameter subset O²y and ox interior), together with the eigenvalue and boundedness constraints therein. Thus SRUW-5-6 are subsumed by Assumption 8.
Hence, for the purposes of Lemma 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-consoldaaa those aligned with the SRUW framework and those introduced carlier for penalized likelihood analysis-into a unified assumption block.
Assumption 13 (Standing Assumptions for proving Theorem 3). The following holds:
1. Identifiability & Smoothness (Assumption 7): the GMM density fetuslyy; cd) is identifiable and E(y;; 0:) = In Jetue (y; ox) is three times continuously differentiable an open ball Bio" ro)-
49

---

2. Compactness & Bounds (Assumption 8): a' is an interior point of compact 0i C [B(,,, ro): mixture weights, means, and covariance/precision eigenvalues satisfy the stated uniform bounds (with Tmin, '), min: 0 o max: min, maax).
3. Data-generating mechanism (Assumption 9): y1,..., yN i..dd from fetuol x²) with latent Za Mult( πT' and yn Zn. k-- N 2
4. Posterior responsibilities (Assumption 10): tu (y; for all Øyv
5. Fisher RE on a cone (Assumption 11); restricted eigenvalue condition for I(x") on C(ca, So) with constant KT >
6. Uniform Hessian concentration (Assumption 12): for radius dn z / ln and
N 2 So In da.
so Ind
sup @e nn AcC)co.S)); || |
sup
A
(N((ā) - A S CH
N
with probability least cf de
7. Penalty penalized coordinates. On the penalized index P, choose X3 that
>S Ao | Vex(((a)))) pe for Ag 1,
egg. a common level Aehosen AoC, V Inda nates excluded from P and from all ||- -lpen, norms.
N with from Lemma Unpenalized coordi-
Assumption (Working high-probability event for ranking analysis). Let Eyma ( lloccpen I Cg /nndd//// be the event in Lemmn 2, and Essc the event on which Equation (21) holds with curvature Kr, 0 over {A C(co, Sa) || S ōR). Under As- sumption 13 and for N large enough, both events hold with probubility at least - SN. where
ōN C d= Crra
+ + di In the proof of Lemma 5, we condition on Eguas O ERSC-
Definition 4 (Quantities for Ranking Consistency). Let a² be the true GMM parameters, and &(X) be the estimator for a given regularization level (with either tied to or fixed).
Noise level.
In da Cg
Anokse
N
the uniform bound on the score at a' from Lemma
2.
Effeeciive curvature for means.
H Ex- (y: a") (12))
where the expectation is with respect
y fetum (::
By Assumption 8, Ex (3:
a²)]
T 2 Tmin and (MF) 2 Omin (since v > 0 and
e, Wie 2 Amin (V))), hence
H (xx) Z Tmin Omin >
KKT remainder for mean coordinates.
Ray (x)( xx,",,",,, a²,,) = Oun Ln (ā(A)) Ope, (x(")) - (((s)()) - mkb)))
We assume that, with high probubility and for the regime (ln
|Rg (á(A), ,>,))) & ERKT Exexr_rom(2 (X), where EXKT (X) depends on (x()= - lx (cf. Lemma 3) and on bounds for Hessian/third derivatives (from Assumptions 7-12). In particular, when || is of order v/ so In da/N.
we take
ExxT (A) S Crem Anoise
for some constant Crem 2 0 over the relevant range of X.
50

---

Noise-screening threshold.
AN (1 Crem + EN) Annises
with small margin eN O.
prese--presevvation threshold. For coordinate j and some ko attaining (or meeting) signal condition for |Piojh
Ag(s;\) X) Ho (a") >/ooool - (1 (1+ Anaise (1 (s) Fe (X);
where Eg > is a fixed margin. Choosing so that A3(;; > O ensures the signal at (ko. j) persists (i.e., Pkay (2) * 0) under the KKT inequalities.
Assumption 15 (True Relevant/Irrelevant Variables). Let Sj denote the set of indices such that at least one component mean has a nonzero j-th entry, i..,, Sp. iff 3k E (1...., .......K) K} with Mkj + O. For j & Sp, we have Pig = 0 for all k. Assumption 16 (Minimum Signal Strength). For each j E si. there exists ko E (l,...,Ky K} such that, for all À up to some detection level A _uppers
Holt (a²) | 2 A (1 6g)
(1+ 11++5) F (A);
(30)
with > 0 fixed. In particular, Equation (30) implies A3(;;A) > X. Moreover, we assume the separation
min JESL.
inf(X 2 NN: where = 0: x Ng(;;\)).
>
(31)
ie., the smallest regularization at which the j-th signal would vanish strictly exceeds the noise- screening threshold Ayy-
Assumption 17 (Regularization Grid). The grid GΛ = ; {Ma)) covers neighborhood of the noise threshold Ayy and extends up to
min X1 with XJ as in Equation (31). JES
Assume p is either tied to X (e.g. x A) or fixed so that its effect is absorbed by constants in the bounds.
Remark 5. The remainder Rky(ā a² in Definition 4 captures the deviation of the mean- coordinate KKT equation from its linearized form. aggregates: (i) off-diagonal Fisher blocks acting on other coordinates of A, (ii) empirical empirical-population curvature fluctuations (Hx I), and (iii) population curvature drift I(ā) I(xx'( from evaluating tions
Under the standing assump- Assumption 12) and the
Assumptions 7-8, Hessian concentration:
Lz-error bound from Lemma 3, one obtains the bound
R&g (²")) , I EKкТ EekT_rem(A) I Crem
on the high-probability event Egnad O ERSC, provided the sparsätysample-size regime ensures Ch V 5o In N= + T.]l@((A)= a=112 1l2 S (e.g., 80 S /N/Inda A crude but sufficient bound for individual coordinates is || | lxx < Al O(v so In da/N) by Lemma 3. Assumption 16 then ensures the effective signal Hots ( (a*)) dominates both the regularization level and the stachastic/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 Egrad D Egsc from Assumption 14, let the grid GX contain points distributed across [0, Àgrid_max max] with
maxx
JES 31
min Xt
x > 0: X -A30;\))
Then, with probability at least
Prank 2
Ma. D -8. where óN == Caraa do
ctt do
MG = 19x,, and D denotes the number of penalized mean coordinates per variable (e.g., K), the following hold:
51

---

1. Relevant Variables. For each j E Si.
Ox(i) 2 Wiigna()) = 4x X A3(:
Let yj == jisj Naignat By the separation Equation (31),
NR > #1XE A An}.
2. Irrelevant Variables. For each j $ Sp
Ok(j) S Nnoise #(XE 9x ANA
Consequently, QR > Nuoises yielding a strict sepuration in ranking scores between relevant and irrelevant variables with probability at least Prank-
Proof of Lemma 5. Work on the high-probability event
Ezank Egrad O ERSC
(on which the uniform gradient bound |ررك1 S Anoise and the KKT-remainder bound (ā(A), at, x)]= < ExKr_Jemi (A) hold simultaneously for all penalized mean coordinates (k,j) and all X E Gx). By Assumptions 13-17 and the union bound,
P(Erank) 2
MM D SN, ón = Ci
de --2 H Hd
+
The ranking score for variable j is
Okli) Σ
I (L...., K} (as(A) + 0)
NEG
KKT expansion. For a penalized mean coordinate (k, j), the KKT condition reads
(A)))
(ay(A) * 0, (k5(A) = 0.
Vu, (v Eujs
Eky E 1,,,
first-order expansion at ax yields
V Vauy IN (ā(A)) V May [x( + Hall (a (T>)( (ks)) + Re)(a(2),
hence
EEg Sks
Hell (oe (T) ()))) REs Raj()(A),
(32)
Relevant variables (j E ST. j E Sy. By Assumption 16 there exists ko with Mias * satisfying the signal condition. Assume, for contradiction, that at a given À E GA we have P (A) = Then ExOJ
[-1, and finas (A) Miaj -Kooj Applying Equation (32) and taking absolute values,
Hott (ax) pien Stoj -
R ((
(a") |MEoo) - | - |Rosl
On Erank
A >
Hell (a") |pãozl
Anoise EKKT (A).
koj
Therefore, whenever
A Helf (aa((")) (a") |pios!
Anoise
EKKT EkKT,rem())
we must have / (A) stricter threshold,
o O. Since A (j: A) 1100,(0" (at) |MGgg) ((+te)))m (1+Eg) Exnes) (A)
< Ag(;AA = (hos(A) #
Hence the count of grid points for which variable j is (at least in one component) active satisfies
Ox ()) Njignali) # GA A A3(;:A)}.
52

---

Taking the minimum over j E Sj gives RR = min Niggnal (7).
Irrelevant variables (j $ (i.). Here Mkj 0 for all k. For given X, the zero solution Pay (2) is KKT-feasible iff
P LN (ā(())) - A.
Using the expansion with nū, = and Pks (A) = 0.
Sk Ray (ā(2), A)] 1 A.
On Syank this is guaranteed whenever
2 Anoise + _KKT____ semm (A).
By definition of the noise threshold Al- = =1++Ew Amosse + (1 (N+ SA EKKT EkKT_rem{ (X') and the monotone/worst-cas domination in its definition, any À > A satisfies X 2 Anoise + EKKT rem (A). Therefore, for each j & Sj. and all À > Aj. all coordinates (An equal zero, so the indicator I(ak: Pag (A) * 0) = 0. Consequently,
Ok(i) S Nnoise == #(A E A AWJ.
Separation and probability. By the separation Assumption Equation (31),
JR
N
Nnoise
yielding strict gap between the scores of relevant and irrelevant variables Eyank- Finally, by the union bound over at most Mg grid points and D) penalized mean coordinates per variable, we obtain
P the above conclusions hold for all 5) 2 1 - -M D D8.
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 and Lemma 5, logether with Theorem for the final KK,,,,, ) choice, hold. 8q = [So] and MNR = p 8o be the number of non-So variables. Moreover, suppose the following assumptions hold:
(a) False negative for So. For any true relevant variable j So and any intermediate set Seur C So with j & Sear the probability of incorrectly rejecting j from S by the BIC criterion is uniformly bounded:
sup P(BICag(i 0eaN) 50) I PN(N).
S... C80
where so pspss) S (N) and (N)
(b) False positive for non-So. For any true non-relevant variable j & So (i.e., j E Un Wo).
given that So has been correctly identified (i.e., Saur So). the probability of BICagli Sn) > 0 is bounded:
P P(BCCag( So) 0) PNNN),
where (N + a N + J0. For BIC, typically (N) = O(N-"s) for some > 0 (e.g., 7n 2 Avmin/ where Avmin 2 is the minimum parumeter-penally gap).
(c) Consistency of BIC-penalized regressions for Rij
So].:
For E Wo. TP ( So] 0) 2 Pogg Pmg((N) (N).
For j E Wo. P( So] Ro()) + 0) 2 1 - Preg (N).
where (wo 2o) Pres (N) ON (1).
53

---

Let the stopping- rule parameter for selecting 3 (and W) be understood so that once S.e. = So holds, each & Sn tested at most times by the BICagf Sn) screening before termination.
For a desired tolerance ES,FP > 0 on the probability of inchuding any false positive in S, it suffices to choose that
c - (l..Fp//MNN) In (1-Px(N))
MNR (1 (1-11-NN(T)) & CS.FP
(33)
For small PN (N), this is well-approximated by
CS.FP
c
MNRPN(N)
In particular, pN decays polynomially in (e.g. N with 7B 2 1/2) and MNRPN(N)
O, then any fixed Cfoed
(e.g., Cfuee or 3) ensures
Pa Plany false positive in S) MNN (1 (1- (1-(1-Px(N)) (N))))) S MNR Cloed PN (NN)N 0.
Then, under (a)-(c) and with chosen according to Equation (33) (or with Gfxed such that MRR Gaxed PN (N) - O), the two-step SRUW procedure recovers the true model structure (Ko, mo, lo. Vo) with probability P(SSccees)) + as N + cc.
Proof. We show each stage succeeds with high probability and then combine them.
Conditioning on ranking accuracy. By Lemma 5, with probability at least - Praitrank lotals the ranking event &aok holds: all so relevant variables in So appear before any pure-noise Wo variables, and Uo appear after Sp. We condition on Eank for the remainder of the argument.
No false negatives in C Fix j E Sn and any Scur C Sn not containing j. By Assumption (a),
P (BBCani( Bear) S 0) I ps(NN.
By union bound over the so true variables,
P P(s g Erank) I so sops(N) =: cs.PN(N) + 0.
Let Es-1OPN be the event So C S. Then PEES-noFN Eyank) 2 1- - ES,.F((N).
No false positives S. On Eank O Es morN we have Seur = Sp. For any j & So, Assumption (b) gives
IP PBICCan So) 0) S PN(N) =:
Let CoF be the event that non-relevant variable is ever added to & after Sp is reached. (A) Distribution-free bound. For any fixed set I candidates examined of any size,
P(Bj J: BICait() So) 0) |J4g
by union bound. In particular,
P PEE Erank n ES.N
< M
which tends to 0 MNR (NN) - 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 consecutive rejections, and that the e decisions in such run (asymptotically) independent more generally, satisfy
54

---

P (n.BICarr So) s o0} 2 (1 - g)
for any e distinct non-relevant candidates (e)) Then the probability that given c-block fails (i.c., contains at least one FP) is
P(block error) Σ 1 g)"
Partition the MN non-relevant indices into [MNn/c] disjoint blocks of size (discarding remain- der needed). By union bound over blocks,
MNR (1 (1- 101-997 918) s MNn (1-(1- .
IP nEES
Erank
O Es EO
S
Therefore, enforcing
In (1-Es.Pp/MNn)
MNR (1-11-NNNV)) CS,PP
c In pal1-pNle (1-Px(N))
yields PEE nop Erank ES-nOFN) I CS,FP. For small PN 1 1- (1 (N)())(N)) cpN (N), so S Es.FP/ (MNPP/MNNPPN
Combining (A)-(B): whenever the independencc/dcccupling condition for runs holds, may the sharper (1 - (1 - 9)" design in the theorem statement; otherwise, the distribution guarantee
MNR PN (N) is valid (and implies the sharper one whenever c is fixed and MN - 0).
Put together $ Sg. Let Eg Es-noHN O Es-noti. Then, with either (A) or (B).
P = So Erank) 2 1 - (e, cs,Fr(N) + 1,
for <s.rr(N) = MNn PN (N) in (A), or (s,Fr(N) MNn (1 -1-1p (o-(NN(N))) in (B).
Consistency of W, Ū, R. Given So, the reverse scan uses BIC-penalized regressions to decide W and R. By Assumption (c),
P W Wo Erak nEG) 2 1- 1-(wo + uo) Preg (N) <w.FP(N) - 1,
and
P By Erank n Eg n W Wo)) > uo Pieg (N) 1.
Final SRUW choice via BIC. By Theorem 2, the final selection of (K, r, E) is consistent, i.e., (PPhina__ooiccc +
Multiplying the stage-wise success probabilities,
P(Overall Success) 2 P(Erank) IP(E9) Erank)
P(W = Wn ...)
P(R
1R ....) choice) 1.
The slowest decaying term among Praitrank Cg (N). (s.F,(N), (wo+ Preg (N), uoPreg (N) governs the overall rate. In particular, under the sharper design (B) with chosen above,
6s,FP(N) MNR (1 (1- (I-PN(N)))) < ES.FP.
55

---

while the distribution-fie design (A) yields csPpp(N) = MNRPN(N) 7 0 whenever MNR PN (N) O.
Theorem (Equivalence of MINAR-SSRUW and MAR on Augmented Data for SRUW). Consider an observation (yn: Ca) where complete-datt data likelihood for Yn given cluster Znk under the SRUW model (K, m, m,rr, r, V) be:
(yn y Un and Cn is its missingness pattern. Let the
yn 2nk = 1;0x) = fu(u an) See(Un Vn R. rreg) findep ( ;imdee)-
Assume an MNARE-SRUW mechunism for missing data, where the probability of the missingness paltern Co, depends only on the cluster membership Enk:
f(Cn Yn, 2nck = 1; 1;²xx) = f(en Znn 1;;b))
II
Eaud
(1 (kd) I-cnd
(34)
Pkd
dESUUUW
where Pkd E (0, 1) are components of the missingness parameter th (thus f(en 2nk = 1; tbk) does not depend on yo)). The observed data under this MNAR2-SRUW model is (ya, Cn), and its likelihood is:
LAMARE-SRUWWW cn; 0))
Wh (yn Zunk = 1;0ks S(co Znk = 1;bb)) dy (35)
Now consider the augmented observed vector Uh (y Cn) Assume Tn is i.i.d. from a mixture model in which the conditional law of the missing part yn is MAR with respect to Yh (i.e., the missingness mechanism does not depend on yn given C,)). Define the **aagmented observed- d likelihood as
(J vn Fuk = 1;00) dyy Me, | Enk = 1; (36)
SMAR 0, (b)
Σ Tk km
where fSRUW(- Znk = 1u) is the same component density as above.
Then, for fixed parameters (0, v²), the observed-date likelihood under MNAR:-SRUW for (y C.) is identical to the likelihood of the augmented observation Th: = (yh. Cn) under the MAR interpretation:
LMNAR--RUW yh Cni O, v/) fMAR 0,
Proof of Theorem 5. Starting from Equation (35), by the MNARz assumption f(cn, Yns Enk
thx) [Cn Znk 1; is independent of y hence
(M C.; 0,bb)
Th (SRUW(y,, "h Enk = =1; f(en Inn = 1; dy
Th f(n 2nk 1;3bk) (ʃ fSSUW(U V" ! 2nzk 1;) dy
Σ Th fu Ink = 1; f(en 2nk 1;bb))
where we set fv( Equation (36) yields
Znk = 1;0g) = - fsRUw( 3.th Znak = 1;;k) dy.. Comparing with
LANAR cai cn; 0, SMAR(Th: 0,4)),
as claimed.
By Theorem 5, the observed-data likelihood under MNAR>-SRUW equals the augmented observed- data likelihood under MAR for every (0, bb) and each observation; hence, the full-sample likelihoods coincide pointwise in (0, v)) The following corollary shows the resulting estimator-level equivalence.
56

---

Corollary (Estimator equivalence under augmentation). Under the conditions of Theorem S. for any fixed (0, 4b and for each observation.
LANARR-SUUW Uh c,; v)) fMA 0, 2/).
Hence the full-sample observed-dee likelihoods coincide, and therefore the sets of maximum likeli- hood estimators for (0, +b) under MNARE-SRUW based on (yy" c) and under MAR the augmented data Uo' = (y²; c) are identical (up to label switching). Moreover, if the same priors (0, ab) are used both formulations, the resulting Bayesian posteriors coincide.
Proof of Corollary 1. By Theorem 5, for every observation and every fixed (0, x/),
LMNAR (Yn Cni 0, xb) SMAR ( 0, tb),
yn (y Cn).
Let D {ly C. ) n denote the sample, assumed i.i.d. under the mixture model as in the theorem. The full-sample observed-data likelihoods are then
N
JANARR.SRIW(DP = II /ANAA Cni 0, ub)
n=1
N
iMAR (D;0,
II n=1
fAR (Th: O,)
By pointwise equality for each factor, obtain the full-sample equality
IMNAR (D; 0 x/) EMAR (D; 0, zb) for all (0, (6).
Consequently, the sets of maximum likelihood estimators coincide:
arg max TNNARL-SRUW 0, = arg max MAR (D; 1b); (0.b)
(9,*)
up to the usual permutations of mixture component labels (label switching), since both likelihoods are invariant under relabelings.
For the Bayesian statement, let II be common prior on (0, v)) admitting density x(0, 2/) with respect common dominating measure. The posteriors are
TNNAR (0, bb D) o (T, 2b) LMNARZ-SRUW(PP vb)
7AAR (0, 2/ ()) (0, ab) LMAR(D) 0 t/).
Since the likelihoods coincide pointwise in (0, v)), the unnormalized posteriors agree, hence their normalizing constants (integrals over the same parameter space) are also equal. Therefore
TMNARZ
D) mAAR (0, 2b D),
again, modulo label switching.
Remark 6. In the MAR interpretation of the augmented data T the components yh have density 1; 0x) in cluster k. and the components Co (which are fully observed within Th have density [(c,, Znk = 1; 1;;b)) in cluster missingness mechanism, given yh and en does not depend 4n itself. The likelihood of ÿh is
fu
The MAR assumption applies y meaning its
formed by integrating out y" from f(yn: Cn. Zukk = 1). which leads to the product fu(ye ) Znk = 1; (x) f(cn]] Znk 1/k)
(. for each cluster k.
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 Y = std((Y) where Y the single-imputed matrix defined in Algorithm which is used for ranking only. In Stage B (role assignment), unpenalized SRUW model under MNARz by exploiting the equivalence to MAR on the augmented data (Y, C). Throughout, V Σ the SRUW partition V (S, R. U, W), and bx denotes EM responsibilities. We write E(-) for (penalized) log- likelikooooso and Q(: 1-) for EM fuuuction...
fit the
57

---

Stage A (ranking): Penalized EM for adaptive-reegulaiized GMM
Penalized objective. likelihood
Given .I...... UN R.,, we maximize the penalized observed-data log-
K
K
(a) Σ log n=
Σ k=1
pA, Ex)
Σ
Pais Maayl (37)
k=1
k=1
Assumption: For each EM fixed (A, the weights Pa computed from the warm start and then held fixed.
At iteration t, with responsibilities
(t) nkk the penalized Q-function is
Qpea(a; (a (1=1)))
log πī.
log Ex
(3h
m4m) To))Tp L mes))
nik
K
|
(38)
1=1
k= 1
E-step.
1)
14g (1-1) 1)(-11
mk
(0)
(4) Thk
Σ
(0)
(39)
Σe= Ti (1-1)
(1r pi (1-1) D 1(1-1))
(t) = n /N. "k
M-step: mixing weights.
M-step: means μ4 with penalty.
Given T ((-1) the subproblem in μr is convex. Let
(8)
(4) nk
m
Then the gradient of the smooth part is
N
lU p4x) TT (1-11)
(z - mm)
(t) 1 ( (+ m (0)).
Tike
=n
n=1
The KKT optimality for coordinate is (t)
n2g
(t- (t-1)
Das
pkj
+
n2k
(2)
"'<4 (wL( (1((1)))
jo Mku
n
(1)
( ((-11)))))) ( (1))s
m
A
coordinate-descent update is
Sx n4
(t) ( (1-1) (1))
(c)
Σ (1-1)) Mko
Mkj
(2) ( (1-1))
Th
ni
)a
where Sx (n) = sugn(u) max { [1u]] X, 03.
M-step: precisions V via weighted graphical lasso. weighted covariance
p (8) fixed, define the responsibility- With
N
S(t))
(gn
2 (0) - (4
(t)
T2k
m=1
Then
arg min --0
leg det v + tr(((() (+) &
v +
2p T
Σ
Pk.ss
M [
with diagonals unpenalized; standard glasso solvers apply.
If P: updated during EM, must add majorization step preserve monotonicity.
58

---

Stage B (role assignment): EM for SRUW under MNARz
Observed likelihood via augmentation.
Let DMAR MNAR [D]. Under MNARZ.
N
( C) =
log
Σ Th M..RMAA( (C.NARrC
a standard mixture on the augmented observation (yh Ca,NNAR).
Complete-data log-likelihood and Q-function.
With θ = (I,, faxx). {zba)) and latent {Enk),
2(0:0
( log EEllog (Y ox.)) E) 2n Sc MNAR (b)))
nk
E-step.
(t-1) (M: (1-1) ((--))((-11)) g((-1)) S (.NNARe MNAR: (
(1)
"
nk
Ee T6
(t-1) R.MARUU: oi (1-1) g(t-
((-11))
M-step: Th and MNARZ parameters.
N
EN (1) EdEDWNAR Cnd 'nk
(1)
(0)
1
(0)
(1)
(1) log
"
N
nk
Pk
Σm (DMNARL
1 (1)
M-step: SRUW data model updates.
Write
fu Org. {) = Setust p
E ( y² a, 3, ()) Sindop " ;7, I).
Let E&[ 5%] denote Gaussian conditional expectations under component k: Then the sufficient statistics are the mixture- moments
E
nk
E&l I- yn Znk =1].
k=1
Cluster block S: for each k, compute
N
N
Fs
(A))))
nk
E L
v.l.
$8.(9
( nk
E&
v.l.
T
n=1
T
n=1
and set 1*
(.(()
Us
(.(()
and 2.11) 2.(1)
k =S. S.100)
yv
5.(1)
(gk
Regression block U let X, = [11 (y Form the mixture--eighteed global moments
N
K
N
K
A
E& X x, ynl.
B
( mk
E
X
y
nk
R=1k 1=
k=
(t) ((
then
A B, and
N
(4)
N
nk
Es
(y
1 - (0)
31(0)) (2 (0) - yn
Y=
1=
Independent block WV:
7(0)
E 19
(
E 1 (2. T 9all 20 (0)).
mk
nk
59

---

Numerical stability. In Stage A, we use warm starts and pathwise (A, p); diagonals unpenalized and enforce 1 > In Stage B, if some ni (f) components per BIC.
tiny, we add a small ridge to Σ or merge/diseard
F
Additional Experiments
F.
Metric Details
We outline some metries that we use to evaluate clustering accuracy and imputation error, which summarized below:
Table 3: Summary of evaluation metrics for clustering and imputation quality
Name
Formula
Range
Imputation error
E (y: gc)2
NRMSE
[0, 1]
(NRMSE, W.)
WNRMSE
[0, 1]
We
Clustering similarity
RI [[R]] max RI) - E[R]]
ARI
[0, 1]
Composite
CIE
a x (1 NRMSE) + / Similarity Score [0. 1]
F.2 Assessing Integrated MNARz Handling
We conduct additional experiment to showcase the performance of handling MNAR pattern within the variable selection framework which the setups are similar [62]. Datasets were simulated with 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 Ynd Σ ô + end, where end ~ A0 1). and ô defined cluster-specific mean shifts with signal strength 2.31 with 811 ô14 822 825 633 δan 7. others zero). For a specific MNAR scenario, the class-specific intercept component v/. was O, and variable-specific slopes b were (1.45, 0.2, 3, 1.45, 0.2, 33. resulting P(M 1|, Za 1) (O))Y)))
In Figure 4 and 5, we plot the boxplots of ARI and NRMSE over 20 replications. We observe that MNAR-bede Ra-asedd approaches attain better results 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, framework demonstrates even slightly higher ARI and lower NRMSE some cases compared to standard MNARz.
F.3 Sensitivity on the choice of
As discussed in [13] and further investigated through Theorem 4, the choice of the hyperparameter plays crucial role in the stepwise construction of the relevant variable set 3, 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 e should ideally be determined by considering the probability of a single incorrect inclusion (px(N), the number of non-relevant variables (MNR), and desired tolerance for overall false positives in & (es,Fr). The experimental results presented here provide practical insights into this interplay and generally support our theoretical conclusions.
60

---

MINARy
MINARZ
MNARYz
1.00
0.50
M Model
0.25
Mean
Mice
missR anger
MNARZ
selvarMNARz
0.75
true mecha
0.50
0.25
0.00
Figure 4: Boxplot of the ARI obtained over 20 replications of simulated data. The theoretical ARIs are represented by red dashed line.
As illustrated Figures and 7, the two datasets (Dataset 14, 80 = 2) show different optimal ranges for This aligns with the theoretical expectation that not a universal constant but rather interacts with dataset-specific factors, which are encapsulated (achieving 100% correct
7,50 = 3; Dataset 2:
by (NN and M. For Dataset 1, the excellent performance with selection of relevant variables) points to very low pv PN(N). This suggests that after the three true relevant variables are found, the remaining four non-Sp variables consistently produce BICai S making small stopping value like 2 both effective and safe. The continued strong performance c 7 (90% correct) further indicates that (NV) small enough that even scanning all MNR non-relevant variables rarely leads to false inclusion. This situation is consistent with theoretical scenario where pv (N) < 1/MNR minimizing the risk of false positives regardless of within reasonable range.
In contrast, Dataset 2, which has non-S, variables (MNR 12), demonstrates clearer trade-
off. Good performance is observed for e = 2 (80%) and = 3 (90%), implying that px (N) is still relatively small. However, the notable decline in performance when So selection) empirically validates a key theoretical concern: if e is set too high and pN is not sufficiently close zero, the algorithm examines more non-So variables while awaiting consecutive negative BICaif values. Each additional variable examined increases the cumulative chance of Type I error non-Sp variable having BICdiff 0 by chance). If the likelihood of obtaining e correct rejections in a row, (1 Px (N))s. diminishes significantly as more variables are processed, false inclusions become more probable. The poor result for
(only 5% correct
in Dataset suggests its pn (N) value makes it unlikely to achieve seven consecutive correct rejections before false positive occurs among the MyR z In in(MNR/a (/w.e/FEFF))) // (N): for given Cs,PP, higher psv (N) larger MN would generally favor a smaller to maintain that error tolerance, or a large fixed e might lead a poorer effective Es.P.. The finding that = 3 is effective for Dataset is consistent with the heuristic used in prior work [13), suggesting it offers practical compromise when px (N) small but non-negligible, and Myr R is moderate. Furthermore, the relative stability of cluster number selection across different
12 candidates. This aligns with the theoretical relationship
61

---

MNARY
MINARZ
MNARYz
1.3
1.1
M Model
0.9
Mean
NRMSE
Mice
missR anger
selvar
1.3
SelvarMNAR2
Figure 5: Boxplot of the NRMSE obtained over 20 replications of simulated data
values suggests that e's main influence on the selection of variables for as predicted by theory, with more indirect effects on the overall model choice, primarily $ is significantly misidentified.
F.4
Computational Times and Scalability
Complexity Analysis
We count arithmetic operations to absolute constants (big-Oh). Throughout this part we denote: N E (samples),
(variables), E N (mixture components). Let MEM be a uniform upper bound on the per-fit number of EM iterations until convergence (Assumption A2 below). In Stage (ranking), the M-step contains graphical-lasse solve per component with at most Mylaso outer iterations. Covariance inverses and log-determinants are computed by Cholesky factorizations.
Assumptions. We work under the following explicit conditions.
A1. The ground-trutt SRUW partition has Derr := |trel + [Umel < D.
A2, Each EM fit terminates in most MEM iterations. For populatioonreegularized EM, this is justified by established convergence rates: either sublinear convergence MAM Op (1/11// (N/ [31] or geometric convergence MEм = (log(1/e)) to achieve e-accuracy [76]. Our analysis uses the more conservative bounded iteration assumption for clarity.
A3. With probability 1 1-od o(1), the Stage-A ranking lists all Def informative variables before any purely irrelevant variables (proved Theorem 3 of the paper).
A4. Let Cylasso (dl denote the arithmetic cost of graphical- lasso solve of size d In gen- eral, Cylasso (d) + d3). Under connected-componce decomposition with largest
block size 8max (as in [73, 41)), Calasso (d) e glassoo Σe v² θ M glasso d 22 _
Per-iteration costs for d-variate GMM. One EM iteration for K-component Gaussian mixture in dimension d has:
62

---

Frequency Metric
Correct Number of Clusters
Correct R Relevant Variables
2
3
100%
87%
8.3%
90%
Proportion of Correct Selections
80%
75%
60%
40%
20%
0%
5
7
100%
100%
98%
81%
80%
75%
60%
45%
40%
20%
5%
0%
Dataset
Dataset_2
Dataset
Dataset
Dataset
Figure 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
E-step: Responsibilities fak require evaluating Gaussian log densitisss for all (n, k). With per-component Σ and log det ΣΦ fixed within the iteration, each evaluation uses a matrix- vector product and quadratic form, costing (df)). Hence ((.)..
M-step (means, weights): Weighted sums are (N)d).
M-step (covariances): Weighted second moments yield O(NK&)). The per-component matrix factorization/iiverrion is (dd)), hence @(Kd))
Therefore one EM iteration dimension d costs
@NNKd)) + O(K.))
and one EM fit (up to convergence) costs
(d) (an()))=(MMM(NNN))) (N + kdd)).
(40)
The classical SRUW selection starts from all D variables and eliminates at time. At step (j = D, D - 1,..., 2), to remove one variable it evaluates j candidates; each evaluation requires an EM fit in dimension j Using Equation (40), the total cost is
D
Σ j CaM(i - = @(MMM Σ
j=2
j=2
Using the polynomial sums:
D² + +1² O(D4)
and
D(D 1)D+ +3D
3D - 1) -(D*)
30
1=1
63

---

Metric
ARI
CIIE
WNRMSE
1.00
0.75
0.50
Average Metric Value
0.25
0.00
1.00
0.75
0.50
0.25
0.00
Dataset 1
Dataset
Dataset 1
Dataset 2
Dataset
Figure 7: Clustering performance of simulated dataset in Section 5 under varying c. We the experiment over 50 replications
we obtain the tight bound
Cteppise (MEM (NKD* + KD*)).
(41)
The D term renders this approach impractical beyond moderate D.
Stage A (Ranking). For each (A, p) on a grid of size Mgrid. we run a penalized EM. Per iteration: the E-step remains (²)) the M-step adds K graphical-lasso solves. By Assumption A4,
per iter cost @(NKD²) O(K Cglasso (D).
Therefore the total Stage-A cost
Crank e MMgia MEM (NKD² + K Calasso (D)))
(42)
Two explicit regimes follow immediately from A4: Cglasso () = + glasso D3)
(Dense/general)
Crank 2 M MEM (NKD² - K K Mglaso Di)):
(43)
(Connected components of size )max) Cglasso (D) = (Mglassc
gglasso Ds
Crank θ Mgris
Mgn MEM (NKD² K MM Ds
(44)
Stage B (Role assignment). A single forward/backward pass evaluates constant number of SRUW- MNARZ fits per newly considered variable. Under A3 the pass stops after Den additions; hence the Stage-B complexity is
Crote = e MEM (NDD2 + KD&))
(45)
Since Den & D (A1), Crste is strictly lower order than Crank and Stage dominates. Combining Equation (43)-Equation (44) and Equation (45), the two-stage complexity is
Cywo-stage Crank + Crole = @ Mana MaM (NKDD K Calasso (D)) + o(Craak).
(46)
64

---

From Equation (41) and Equation (46), the speedup factor satisfies
Cstepwise
NKD + KD
Ctwo-stage
Marid (NKD2 + C Classoo (D)
Two explicit lower bounds follow.
If NTNE 2 Mglasso (E-step dominates Stage A): using Equation (43),
Cstepwise Ω
D2
Cwo-stage
Mgrid
If N < Mg Mglasoo D (glasso dominates Stage A): using Equation (43),
Cstepwise
D2
Ω
Cwo-sttge
Maria Mglusso
Under comnected-componett sparsity with block size Smax Equation (44) the second case strengthens
to
Cstepwise
D2
N
})
Ω
min 1,
Cewo-stage
Mard
Mglasso maax
In all regimes the speedup is 12 / MMgna), and strictly larger when sparsity (small Smax) is present.
Empirical Validation
We report wall-clock times (seconds) for representative scenarios, confirming the predicted polyno- mial speedups:
Scenario
N
D
K
SelvarMNARz (s)
Clustvarsel (s)
Speedup
Varying D
750
15
12.2
190
~55
Varying D
750
21
4
14,4
640
444 x
Varying D
750
27
15.5
2054
132x
Varying N
1000
20
4
19.0
838
444
Varying K
750
20
12
9.77
1242
127x
These measurements consistent with the theory: the two-stage method scales like D in the dense case (and better under sparsity), whereas backward stepwise scales as D
Dependency Computational Times under Degree of Missingness
We analyze the computational dependency of our framework on the proportion of missing data, revealing an advantggcous 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 X,, denote the number of observed entries in observation y.. The E-step cost for Gaussian mixture models scales as:
N
K
Cr.step
Σ
X2
=0 =K(NE(X2)
=
1=1
For different missingness mechanisms:
MCAR at rate r: X ~ Binom D, r), yielding
B|x Var(X) + E[X] D(1-r)r+ D² r)2
This produces quadratic decrease in E-step work as 1.
65

---

The same complexity bound applies. with X,, representing the random count of observed coordinates per observation. Each per-record Gaussian density evaluation scales quadratically with observed entries, maintaining the e (KKNE[X complexity.
MAR/MNARZ:
Stage A employs efficient single imputation and remains largely insensitive to missing rates, while Stage inherits the beneficial E[X²) scaling. Consequently, total runtime decreases monotonically with increasing missingness rates.
We empirically verified this property by fixing D2) (1000, 14, and systematically increas- ing the missing data rate under two mixed missingness scenarios. The results, summarized in Table 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
Algorithm Runtime vs. Degree of Missingness
60
Runtime (seconde)
30
40%
Per-Variable Missing Data Rate
Missingness Mechanism Mixed (MAR-MNAR)
Mixed IMCAR-MAR--NARR
Figure 8: Runtime the proportion of missing data increases.
Three key observations emerge:
Runtime consistently decreases with higher missing rates, as predicted by the B[X"] scaling. The MNARz block processes smaller observed patterns during E-step calculations, reducing computational burden for conditional expectations and complete-data sufficient statistics.
2. The slight variation at 30%-50% missingness in the ftree-mechanish scenario likely stems from random allocation of missing positions. When missing data concentrates in clustering variables (S), the number of EM iterations for convergence may vary, the overall inverse relationship persists.
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.
66

---

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 dataset of D 14 variables, where the true clustering variable set is S' (1,2) 2) (remaining variables play roles of TR, U, W according to scenario in section We generated missing data using mixture of mechanisms:
MCAR: Some variables have missing values completely at random
MAR: Missingness depends only on observed components yn
MNARy: Missingness depends directly on unobserved values yM (actively violating MNARZ assumption)
The MNARy mechanism (see [62] for precise formulation) specifically tests our model's robustness to misspecification, it depends 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. MNARz + SelvarMix: Estimate GMM-MNARZ model and impute missing data, then run SelvarMix [13] on the imputed dataset
2. Multiple Imputation variants: Using both missRanger and geimputeR (cite here) with Mclust
3. VarSelLCM: competing joint modeling approach
Table Experiment with mixed MAR+MNARY (True set {1,2])
Method
ARI
NRMSE
Relevant Variables Selected
SelvarKNARz (ours)
0.808
0.176
1, 2 (correct)
VarSelLCM
0.779
1-11 (extra variables)
missRanger Mclust
0.799
0.190
geimputeR Mclust
0.774
0.413
MMARZ SelvarMix
0.328
0.070
(missing variable)
Table Experiment with mixed MCAR+MAR-MNARS (True set (1,2))
Method
ARI
NRMSE
Relevant Variables Selected
SelvarKNARz (ours)
0.808
0.181
1, 2 (correct)
VarSelLCM
0.772
1-11 (exra variables)
missRanger Mcluat
0.783
0.198
geimputeR Mclust
0.755
0.401
MMARZ SelvarMix 0.328
0.087
(missing variable)
The poor performance of the MNARZ SelvarMix baseline prompted 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 Mixed (MAR+MNARy) Data (Rmixmod backend)
Method
ART
NRMSE
Relevant Variables
SelvarMNARz (Ours)
0.778
0.162
1,2
Decoupled (SI MNARz)
0.328
0,263
Decoupled (KI mi missRanger)
0.336
0.257
Decoupled (HI MNARz)
0.328
0.331
67

---

Table Performance on Mixed (MCAR-MAR+MNNRR)) Data (Rmixmod backend)
Method
ARI
NRMSE
Relevant Variables
SelvarMNARz (Ours)
0.773
0.164
1.2
Decoupled (SI MNARz)
0.328
0.387
Decoupled (MI MNARz)
0.336
0.322
Decoupled (MI missRanger)
0.328
0.248
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:
Uncertainty Propagation. Decoupled pipelines treat completed data observed, discarding posterior uncertainty in y Our EM-based framework propagates this uncertainty through all parameter and role updates.
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 empirically more stable 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 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 hierarchical clustering (HC) based initialization (Ward's linkage on Euclidean distances with cluster centers extracted from the dendrogram cut K) against multiple random initializations (MIs). HC places initial centers in high-density regions, yielding more stable responsibilities at the first E 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 (higher ARI is better). MIs: Multiple random initializations.
Scenario
MIs (ARI)
HC (ARI)
I
0.283
0.317
2
0.479
0,533
3
0,551
0.536
0.441
0.654
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.
F66 More on Transcriptome Dataset
Background and Prior Analyses
Dataset. We analyze 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], retain all 1267 genes: 1149 are fully observed: 118
68

---

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.. cluster of 15 genes co-clustered with 4 well-studied markers). SelvarClustVV [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-relese relevant, together with PI-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. 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 E .,....2, with = 5, spectral distance weights Pk. and P& LC structure, the selected model yields 18 clusters and a global role assignment with P1-P4 in S and P5-P7 in U. This agrees with prior work on P5 (not clustering-relevan) relevant) but differs by reclassifying P6-P7 from $ (in [37, 39J) to U.
Cluster-level diagnostics clarify the difference. Let R2 denote the coefficient of determination from regressing the U-block (here, P5-P7) on the S-block (P1-P4) within cluster Table shows pronounced heterogeneity:
Several clusters (e.g., 6, 7, 8, 10, 12, 18) have R2 > 0.60, indicating that P5-P7 are largely explained by P1-P4 within those clusters. For these groups, assigning P5-P7 to U is appropriate.
A large aggregate (Cluster I) exhibits low R² alongside flat S-profiles and detectable UJ- activity, suggesting local decoupling of P5-P7 from P1-P4. This mirrors the biological intuition behind earlier inclusion of P6-P7 in S.
Hence, both the earlier global S-assignment (P6-P7 in S) and our global U-assignment (P6-P7 in U) 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-specifie,
Why does our global BIC prefer P5-P7 in U?? Two factors are play:
1. Global model selection. BIC aggregates fit-complexity tradeoff's across all clusters. Because many clusters exhibit high R (P5-P7 explained by P1-P4), the global criterion prefers a parsimonious $ (P1-P4) and assigns PS-P7 to U.
2. MNARZ identifiability and shrinkage. Under MNARZ, class-specific missingness param- eters Pkd absorb class-linked absence patterns. When a project effectively behaves as MAR/MCAR for many clusters, the fitted Pkd across k becomes nearly homogencous and the conditional dependence of P5-P7 on P1-P4 becomes tighter, further favoring U globally.
Biological reading consistent with both views. P2 (iron signaling) is reaffirmed core driver. P5-P7 behave as late/stress outputs strongly coupled to P1-P4 in many clusters (high RC), but retain independent variation in sizable subset (e.g., Cluster 1). This reconciles the prior decision place P6-P7 in S (emphasizing sypocotyl-centric signals) with our global U assignment (emphasizing aggregate parsimony across all clusters).
Practical implication. Our diagnostics suggest a natural extension: cluster-adaptive role assignment or hierarchical prior tying per-cluster roles, which would allow P6-P7 to enter S only where local evidence (low residual error) warrants while retaining parsimony elsewhere. This aligns with the heterogeneity revealed by R² and preserves the strengths of both global viewpoints.
Our unified S/U 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-bar R-baedd analyses.
69

---

Terk SEAAET
Cleorr
Ckatter
Clptee
Clucore
Chackers
Clacte
Cluoe7?
Chator
Elaster
Clacter
Cliter1
Claiter
Clacoe
Clacter
Figure 9: Mean expression profiles for each of the 18 clusters, displayed in separate panels. Light region indicates irrelevant P.
G Additional Details on Related Work in the Literature
G.1 Model-based Clustering
Model-based clustering conceptualizes clustering 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-bart based methods, this paradigm enables principled inference, allowing for parameter estimation via likelihood- techniques and objective model selection to determine the number of clusters.
prototypical instance the GMM, wherein each cluster is characterized by 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, 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 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 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 address this issue [65]. A notable advantage of the Bayesian approach is the incorporation of priors model complexity, e.g., via Reversible Jump MCMC or birth-death processes, to infer the number of clusters. Despite their ficxibility and robustness, fully Bayesian clustering methods can be computationally demanding, particularly when MCMC chains converge slowly or require extensive post- 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 AIC [1] or BIC [58]. While effective for moderate- dimensional settings, these methods become computationally prohibitive as the number of features D increases.
70

---

The advent of penalized likelihood methods improved scalability and enabled sparse modeling. The LASSO [70]. a seminal technique using an L, penalty, enables both variable selection and coefficient shrinkage. Predecessors include the nonnegative garrote [9] and ridge regression [28], the latter using an 1/2 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, vacc, 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 heuristies in the clustvarsel R package streamline EM evaluations.
Expanding this idea, Maugis et al. [38] proposed 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, 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 SRUW 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 hooeehherrstiicadiussed 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 (MJ) [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 im- plementation, fits univariate models conditionally and iteratively imputes missing values, supporting mixed data types and nonlinear relationships.
High-dimensional settings pose 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 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
71

---

like Fed-MIWAE [4] for and are best suited for large sample sizes, though interpretability remains a challenge.
contexts. These models assume latent manifold structure
Addressing MNAR data in clustering models remains difficult due to identifiability issues and the need for strong assumptions auxiliary information. Two general strategies exist: (1) Selection models, specifying 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 augment clustering models with missingness-inee -nforremss---de informed constraints to improve identifiability.
72


Reviews:
[{"review_id": "sSkrA5kIDZ", "reviewer": "Reviewer_BwuH", "summary": "This paper presents a novel unified framework for model-based clustering that effectively tackles both variable selection and missing data (specifically, MNAR) simultaneously. By integrating a data-driven penalty and explicitly modeling missingness with latent class membership, the proposed method, estimated via an extended EM algorithm, demonstrates theoretical guarantees of consistency and empirical improvements in identifying informative variables and clustering complex datasets.", "questions": "### Mixture of missing type\n\nWhile the proposed framework demonstrates strong performance under the assumed MNARz mechanism, its applicability to real-world datasets often characterized by an unknown mixture of MCAR, MAR, and diverse MNAR patterns could be further clarified. To enhance the paper's comprehensiveness and provide additional insights for practitioners, it would be highly beneficial to include empirical results on synthetic datasets that deliberately incorporate a blend of these various missing data patterns. Such an analysis would offer valuable understanding into the algorithm's robustness and generalizability beyond a single assumed missing data mechanism.\n  \n### Conclusion from Figure 1\nThe authors summarize the first synthetic experiment as \"Its performance declines only modestly as missingness increases, while the impute-then-cluster baselines deteriorate sharply,\" The figure is quite small but it looks like all four methods deteriorate in the same way, yet the proposed algorithm is definitely has more tolerance on missing ratio. It would be great to enlarge this figure and ensure clearer differentiation between the performance curves. Some statistical tests to quantitatively supporting the claim is welcomed to report.\n  \n### Decoupled versus joint approaches\nThe proposed framework admirably integrates variable selection with missing not at random (MNAR) modeling by extending the SRUW framework under the MNARz assumption. However, since the core innovation builds upon the MNARz mechanism, it would significantly strengthen the empirical contribution if the authors compared their method against a more explicit baseline that combines vanilla MNARz modeling with standard variable selection applied to fully observed data. Specifically, I encourage the authors to construct a hybrid baseline that uses:\n\n- MNARz modeling (e.g., a GMM with cluster-dependent missingness probabilities as in [59]) for handling missingness,\n- Followed by variable selection methods (e.g., penalized likelihood via adaptive LASSO or SRUW role assignment) applied to the internally imputed or completed data.\n\nThis would help isolate the gain provided by joint modeling of selection and missingness versus decoupled approaches. Even a simplified version of such a pipeline would provide valuable insight into whether the proposed method’s improvements stem from its structural innovations or primarily from leveraging MNARz assumptions.\n\t\nAsking for comparison with vanilla MNARz + variable selection isolates whether joint modeling (i.e., their contribution) gives significant advantage over a two-step pipeline.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "### Strengths\n\n- The paper integrates two significant challenges in model-based clustering: variable selection and handling MNAR data. This unified approach is a notable advancement, as most existing methods typically address these issues in isolation.\n- The framework provides strong theoretical backing, including proofs for identifiability, asymptotic consistency, and selection consistency, even in the presence of missing data.\n\n### Weaknesses\n\n- As acknowledged by the authors, the current framework is limited to continuous data, specifically via Gaussian mixture models. This restricts its applicability to categorical or mixed-type variables, which are common in many real-world datasets.\n- While powerful for targeted MNARz scenarios, the framework currently does not offer an explicit mechanism for automatically identifying underlying missing data patterns in mixed real-world settings. This aspect might be a practical consideration for its broader application in highly heterogeneous missingness scenarios.\n- In the synthetic experiments the MAR and MNAR missingness regimes are simulated separately — not simultaneously. For practical usage, it would be highly valuable to see empirical results on synthetic datasets that deliberately incorporate a blend of different types or complexities of MNAR patterns. This would offer crucial insights into the algorithm's performance and robustness when facing the nuanced and heterogeneous MNAR scenarios commonly found in real-world applications.", "was_revised": true, "final_justification": "As per the system’s instructions, paste my prior final justification:\n\nThe new results support the view that the joint approach benefits not only from its statistical formulation but also from the stability of its integrated estimation process.\n\nTaken together with the discussion raised by other reviewers, these findings reinforce my original assessment: the paper offers a meaningful contribution, and the authors’ thoughtful responses have further strengthened its empirical and practical value.", "rebuttal": "We are deeply grateful to you for your positive assessment and for providing such insightful and constructive feedback. We are delighted that you recognized the strengths of our unified approach and its theoretical backing. The suggestions for improving the empirical evaluation are excellent, and we agree that incorporating them will significantly strengthen the paper's contribution and practical relevance. We will address your concerns and suggestions below\n\n## **1. On handling a mixture of missing data types**\nOur current experiments focused on \"pure\" MAR and MNAR scenarios to clearly isolate and demonstrate the benefits of our MNARz-aware model. You are correct that a mixed-mechanism scenario is a more realistic and challenging test of robustness. We will first demonstrate the theoretical behavior and our new small simulation study:\n+ Our framework is ultimately designed to handle the most challenging case (MNAR). When MCAR or MAR data are present, our model is technically misspecified, but it is a \"safe\" misspecification. The MNARz model can be seen as a more general class that can approximate MAR/MCAR. For instance, if a variable is MAR, its missingness does not depend on the latent class $z_n$. In this case, the EM algorithm should ideally learn similar missingness probabilities $\\rho_k$ across all clusters $k$ for that variable, effectively recovering a MAR-like mechanism. While not perfectly efficient, the model should not fail catastrophically.\n+ To address this question empirically, we have conducted a small simulation study suggested by the reviewer. We generated a synthetic dataset with 14 variables where $\\{1,2\\}$ is the true relevant set. To account for mixed missingness mechanisms, a subset of variables is made missing via an MCAR mechanism, another subset via MAR, and a final subset via our MNARy mechanism. Note that the MNARy mechanism indicates that the missingness probability depends on the unobserved value $y_{nd}$ itself which is also represent the misspecification of our framework since our model is currently built on MNARz.\n+ We denote \"MNARz + SelvarMix\" to be the suggested baseline from the reviewer. More details are in answer for question 3.\n+ Based on the results, our SelvarMNARz framework still demonstrates superior performance in both clustering accuracy (ARI) and, crucially, in correct variable selection compared to MAR-based methods (like SelvarClustMV or impute-then-cluster pipelines) even when the data contains a mixture of missingness mechanisms (including patterns that misspecify our model's assumption)\n+ The reason is that it can correctly model the most destructive component (the MNAR part), while being robust enough to handle the simpler MAR/MCAR patterns.\n\nHere is the result:\n\n**Table 1: Performance on Mixed (MAR+MNARy) Data**\n\n|Method              |   ARI| NRMSE|Relevant Variables|\n|:-------------------|-----:|-----:|:--------------------------------------------|\n|SelvarMNARz (Ours)  | 0.808| 0.176|1, 2 (Correct)|\n|VarselLCM           | 0.779|    - |1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 (Incorrect)|\n|missRanger + Mclust | 0.799| 0.193|-|\n|gcimputeR + Mclust  | 0.774| 0.413|-|\n|MNARz + SelvarMix   | 0.328| 0.070|1 (Incorrect)|\n\n**Table 2: Performance on Mixed (MCAR + MAR + MNARy) Data**\n\n|Method              |   ARI| NRMSE|Relevant Variables|\n|:-------------------|-----:|-----:|:--------------------------------------------|\n|SelvarMNARz (Ours)  | 0.808| 0.181|1, 2 (Correct)|\n|VarselLCM           | 0.772|     -|1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 (Incorrect)|\n|missRanger + Mclust | 0.783| 0.198|- |\n|gcimputeR + Mclust  | 0.755| 0.401|- |\n|MNARz + SelvarMix   | 0.328| 0.087|1 (Incorrect)|\n\n### **2. On the Clarity and Interpretation of Figure 1**\nWe thank the reviewer for this observation and agree that our original presentation and description were not clear enough. We will list out the clarification below:\n+ Our intended claim was not that other methods do not deteriorate, but that the rate of deterioration is much sharper for the impute-then-cluster baselines, especially under the MNAR mechanism at high missingness rates (30-50%). Our method, while also declining in performance (as is inevitable), maintains a significantly higher level of performance, showcasing its greater tolerance.\n+ You are right that the figure is dense, and we appreciate you pointing this out. In preparing the submission, we were balancing the need for comprehensive results against strict page limits, which led to the current compact layout. Our intention was for the high-resolution vector format to allow for detailed inspection by zooming in on the electronic version. However, we recognize this is not always convenient. We will therefore redesign the figure in the revised manuscript to be immediately legible and ensure all performance curves to be visible.\n+ Based on the reviewer's feedback, we have performed pairwise one-sided Welch's t-tests to compare the performance of our proposed SelvarMNARz framework against all other baselines. We focused on the most challenging scenario discussed in the paper: MNAR data at a high 50% missingness rate. The results, including Bonferroni correction for multiple comparisons and 20  replications (sample size), are summarized in the table below.\n\n**Table 3: Welch's t-test for Mean ARI between SelvarMNARz and Clustvarsel with $\\alpha = 0.05$**\n\n|Model       | Mean ARI| Std. Dev.|Comparison      |p-value |Signif. |\n|:-----------|--------:|---------:|:---------------|:-------|:-------|\n|SelvarMNARz |    0.511|     0.052|-               |-       |NA      |\n|Clustvarsel |    0.363|     0.088|vs. SelvarMNARz |<0.001  |****    |\n|Selvar      |    0.348|     0.108|vs. SelvarMNARz |<0.001  |****    |\n|VarselLCM   |    0.344|     0.101|vs. SelvarMNARz |<0.001  |****    |\n\n+ Based on the p-values, the tests support for our central claim. This demonstrate that the satisfactory clustering performance (as measured by ARI) of our method is not an artifact of experimental randomness but a statistically significant and substantial improvement over all baseline methods.\n\n### **3. On Decoupled versus Joint Approaches**\nWe appreciate for a highly insightful suggestion. We agree that comparing our joint model to a strong, two-step (decoupled) baseline that also uses MNARz is the best way to isolate the benefits of our integrated approach. We have conducted this new experiment as requested. Below we briefly summarize this newbaseline which we called \"MNARz + SelvarMix\":\n\n+ Step 1 (Imputation): We first apply a standard MNARz-GMM model (as in [59]) to the incomplete data. This model is used solely to produce a single, complete imputed dataset based on the final conditional expectations.\n+ Step 2 (Variable Selection): We then apply our variable selection and role-assignment procedure as in Celeux et al., 2019 (SelvarMix) to this newly completed dataset.\n\nThe results for this new baseline are included in the tables in Section 1 above. The findings are stark: the decoupled MNARz + SelvarMix pipeline performs very poorly, with a low ARI and incorrect variable selection. This empirical result can be explained in a more theoretical manner. The decoupled approach suffers from a **fundamental flaw**: it fails to propagate uncertainty from the imputation step to the selection step.\n\n+ The Decoupled Approach: The imputation in Step 1 produces a single \"best guess\" for the missing values. This completed dataset is then treated by the selection algorithm in Step 2 as if it were the true, complete data. This ignores the fact that the imputed values are uncertain estimates. The variable selection model in Step 2 becomes overconfident in the imputed data, leading to biased parameter estimates and incorrect selection of variables. The low NRMSE of this baseline is deceptive; it indicates stable imputations (low variance), but these imputations are based on a misspecified model that ignores the true variable roles (high bias), leading to poor downstream clustering and selection.\n+ Our Joint Approach: Our integrated framework avoids this pitfall. By performing imputation and selection jointly, the uncertainty about the missing values (represented by their posterior distributions) is naturally and correctly propagated through every parameter update. The variable selection process is always aware of the missingness, leading to more robust and accurate parameter estimates and, consequently, superior clustering and variable selection performance."}, {"review_id": "ISVewVektf", "reviewer": "Reviewer_rfVw", "summary": "The paper introduces a variable selection framework in model-based clustering that handles data that is missing not at random (MNAR). The paper extends the SRUW model used in variable selection for model-based clustering in high-dimensional settings using the missing not at random model called MNARz which takes into account which class the data is coming from. They propose an EM style algorithm combined with a transformation of MNARz to MAR. The paper shows that under specific assumptions the proposed framework maximizes identifiability and consistency of the SURW framework and with high probability selects a consistent set of relevant features. They also experimentally show that the methods modified with MNARz handles missing not at random better than the standard methods.", "questions": "1. How dependent is the framework on the MNARz scheme? If the missingness is slightly different from MNARz (maybe for instance if it also relies on some other parameter as well), does that break things? Does it break things in practice?\n\n2. How do you initialize the EM algorithm? Is a random starting point good enough or do you need to prime it in a specific way?\n\n3. How does the computation scale with parameters such as $D,N$ and $K$? How scalable is the algorithm? What is the dependency of the algorithm's complexity on the degree of missingness?", "limitations": "The paper highlights the limitations such as the restriction to Gaussian mixture models in continuous setting, in the discussion.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n\n1. The paper explores the variable selection in model-based clustering. They introduce a method that extends SRUW model to handle missing not at random data using the MNARz method. They show that the method they introduce satisfies identifiability and consistency under the SRUW model and with high probability, can extract a set of reliable and consistent features.\n\n2. The algorithm is generally easy to follow and uses an EM based approach. The overall organization of the paper and the arguments are clear. The paper is overall well-written.\n\n3. The paper shows the practicality of the MNARz based approaches and shows experimentally how they do against standard methods.\n\nWeaknesses\n\n1. It is unclear how the computational complexity of the algorithm scales with the parameters such as $D,N$ and $K$. \n\n2. The conditions for the correctness of the algorithm are mentioned in the supplementary material. It might be better to have at least some informal statement about them in the main body because they would be important in ascertaining the validity of the theorems.", "was_revised": true, "final_justification": "The authors addressed my concerns about scalability, degree of missingness as well as on the dependence of the missingness model and the algorithm initialization.", "rebuttal": "We sincerely thank the reviewer for their positive and encouraging feedback. We are delighted that you found our paper well-written and our framework's contributions to be significant. The questions raised are very insightful, and we will address each question on the robustness, initialization, and scalability of our method in detail below.\n\n### **1. On the Dependence on the MNARz Scheme and Robustness**\nTheoretically, our consistency and identifiability guarantees (Theorems 1-3) hold when data is generated from the specified SRUW-MNARz model. Severe model misspecification would formally invalidate these guarantees, as in most statistical models.\n\nIn practice, however, models are rarely perfectly specified. The key question is whether our framework \"breaks\" gracefully. We have conducted experiments that directly address this. In Appendix E.2 (Figure 4), we tested our method on data that violate the pure MNARz assumption:\n\n1. MNARy: The missingness probability depends on the unobserved value $y_{nd}$ itself.\n\n2. MNARyz: The missingness probability depends on both the unobserved value $y_{nd}$ and the latent class $z_n$.\n\nWe also performed a small study under the case mixture of missing mechanisms as suggested by reviewer BwuH. Our framework remains competitive and consistently outperforms MAR-based baselines, because the latent-class link captured by MNARz explains a major portion of informative missingness. This also aligns with the findings of Sportisse's work, which we use for modelling MNAR data, cited as [59] in our paper.\n\n### **2. On the Initialization of the EM Algorithm**\n\nEM for finite mixtures is famously sensitive to initialization. So, we use a standard and effective multi-step initialization procedure:\n\n1. We perform a simple, fast imputation of the missing data to get a complete dataset. The imputation is used because the ranking stage is essential to avoid misclassification of variable roles.\n2. We apply hierarchical agglomerative clustering or $K$-means on this completed dataset to obtain an initial partition of the data into $K$ clusters.\n3. From this initial partition, we compute the starting values for all our model parameters: mixing proportions ($\\pi_k$), means ($\\mathbf\\mu_k$), and covariance matrices ($\\mathbf\\Sigma_k$). ($\\alpha_k, \\beta_k$) are also initialized. These serve as the $\\mathbf\\theta^{(0)}$ for the main EM algorithm.\n\nThis hierarchical clustering (HC) based initialization provides a much more sensible starting point than random initialization, typically placing the initial cluster centers in high-density regions of the data and greatly improving the stability and convergence of the EM algorithm. Our results (see `Table 4`) confirm it is more robust than multiple random starts.\n\n**Table 4: Comparison of EM Initialization Methods under 8 Scenarios**\n\n*MIs: Multiple Initializations*\n\n|Scenario|Initialization Method|ARI|\n|:-|:-|:-|\n|1|MIs|0.283|\n||HC|0.317|\n|2|MIs|0.479|\n||HC|0.533|\n|3|MIs|0.551|\n||HC|0.536|\n|4|MIs|0.441|\n||HC|0.654|\n|5|MIs|0.760|\n||HC|0.771|\n|6|MIs|0.711|\n||HC|0.778|\n|7|MIs|0.716|\n||HC|0.780|\n|8|MIs|0.785|\n||HC|0.786|\n\n### **3. On Computational Complexity and Scalability**\nWe thank the reviewer for their crucial feedback.\nWe address the valid concerns about scalability by comparing our two-stage, Lasso-like procedure proposed by [2] and adapted in our work, with the classical backward stepwise SRUW version from [3]. First, some terms and key assumptions are reiterated:\n\n1. The true underlying model is sparse. The number of relevant variables, $|\\mathbb{S}\\_\\text{true}|$, and redundant variables, $|\\mathbb{U}\\_\\text{true}|$, define the effective dimension $D\\_\\text{eff} = |\\mathbb{S}\\_\\text{true}| + |\\mathbb{U}\\_\\text{true}|$, which satisfies $D\\_\\text{eff}\\ll D$.\n\n2. The EM algorithms used in both methods converge to a solution within a reasonable number of iterations, denoted as $M_\\text{EM}$. Note that we can assume sublinear convergence (e.g., [1]), typically $T = \\mathcal{O}_p(1/\\sqrt{N})$ or the population EM converges in $T = O(log(1/\\epsilon))$ iterations such that the output are $O(\\epsilon)$-close to the ground truth for spherical Gaussian mixture model. Still, we believe our assumption is more conservative and easier to illustrate the computational advantage.\n\n3. We assume the variable ranking stage is consistent, as supported by the theoretical results in our paper (Lemma 2 in our manuscript). This implies that, with high probability, all $D_\\text{eff}$ informative variables are ranked before the uninformative variables.\n\n4. We define $M_\\text{grid}$ as the grid size for $\\lambda$ used in the ranking stage, and $d$ as the number of dimensions in a model during an intermediate step of an algorithm, where $d < D$. Let $M_\\text{Lasso}$ be the number of outer-loop iterations for the glasso solver to converge.\n\n#### **3.1 Classical Backward Stepwise SRUW**\n\nThis method starts with all $D$ variables and iteratively removes them. At each step with $j$ variables, it must fit and evaluate $j$ different models to decide which variable to remove.\n\n* Fitting a model with $d$ variables costs $\\mathcal{O}(M_{\\text{EM}}\\cdot (NKd^2+Kd^3))$.\n\n* The backward stepwise algorithm uses this model fit as a building block. At each step $j$ (from $j=D$ down to 2), the algorithm must fit and compare $j$ different models, each with $j-1$ variables, to decide which variable to eliminate.\n\n* The total removal cost is:\n$$\nC_{\\text{stepwise}}=\\mathcal{O}(M_{\\text{EM}}\\cdot (NKD^4+KD^5))\n$$\n\nThe $\\mathcal{O}(D^5)$ complexity makes this approach computationally infeasible for even moderately high dimensions (e.g., $D>100$).\n\n#### **3.2 Two-Stage Lasso-Like SRUW**\n\nOur proposed method avoids this combinatorial explosion.\n\n* Stage 1: Variable Ranking: This is the main computational bottleneck. It involves running a penalized EM algorithm over a grid of $M_{\\text{grid}}$ regularization parameters. The cost is dominated by E-step and graphical Lasso updates within the M-step, scaling with $D^2$ (see [4], [5]).\n$$\nC_{\\text{rank}}=\\mathcal{O}(M_{\\text{grid}}\\cdot M_{\\text{EM}}\\cdot(NKD^2+KM_{\\text{glasso}}D^2))\n$$\n\n* Stage 2: Role Assignment: This stage performs forward and backward scans on the ranked list to identify relevant, redundant, and independent variables. These steps operate on a small number of variables ($|\\mathbb{S}_\\text{true}|$, $|\\mathbb{U}_\\text{true}|$, etc.), so its complexity is $\\Omega(D_\\text{eff})$, not $D$.\n\nUnder the sparsity assumption ($D_{\\text{eff}}\\ll D$), the overall complexity is:\n$$C_{\\text{lasso}}\\approx C_{\\text{rank}}=\\mathcal{O}(M_{\\text{grid}}\\cdot M_{\\text{EM}}\\cdot K\\cdot D^2\\cdot (N+M_{\\text{glasso}}))$$\n\n#### **3.3 Computational Speedup**\n\nThe speedup is the ratio $\\frac{C_{\\text{stepwise}}}{C_{\\text{lasso}}}$.\n\n* If $N>D$, the speedup is $\\frac{\\mathcal{O}(NKD^4)}{\\mathcal{O}(NKD^2\\cdot M_{\\text{grid}})}=\\mathcal{O}(\\frac{D^2}{M_{\\text{grid}}})$.\n* If $D>N$, the speedup is $\\frac{\\mathcal{O}(KD^5)}{\\mathcal{O}(NKD^2\\cdot M_{\\text{grid}})}=\\mathcal{O}(\\frac{D^3}{N\\cdot M_{\\text{grid}}})$.\n\nIn all practical scenarios, our two-stage procedure is polynomially faster, with a speedup of at least $\\mathcal{O}(D^2)$. To validate this theoretical prediction, we ran a direct, head-to-head comparison on a focused set of scenarios. The results, combined with previous analysis, demonstrate that our method is a computationally feasible and powerful tool for the high-dimensional problems it is designed to solve.\n\n**Table 5: Empirical Runtime (in seconds) and Speedup Factor**\n\n|Scenario|N|D|K|SelvarMNARz (s)|Clustvarsel (s)|Speedup Factor|\n|:-|:-:|:-:|:-:|:-|:-|:-|\n|Varying D|750|15|4|12.2|190|~15x|\n|Varying D|750|21|4|14.4|640|~44x|\n|Varying D|750|27|4|15.5|2054|~132x|\n|Varying N|1000|20|4|19.0|838|~44x|\n|Varying K|750|20|12|9.77|1242|~127x|\n\n### **4. On the Dependency of Complexity on the Degree of Missingness**\nWe have conducted a new empirical study to provide a definitive answer.\n\nOur two-stage framework's runtime has two components: a ranking stage on imputed data and a role-assignment stage on the original incomplete data, whose cost is theoretically dependent on the average number of observed variables. To test this, we fixed the problem size ($N=1000,D=14,K=4$) and varied the missing rate from 20% to 80% across two complex, mixed-mechanism scenarios. The empirical results, presented in the table below:\n\n**Table 6: Empirical Runtime (in seconds) with Increasing Missing Rate under Mixed Missingness Mechanism**\n\n|Missingness Mechanism|20% (s)|30% (s)|50% (s)|80% (s)|\n|:-|:-|:-|:-|:-|\n|Mixed (MAR+MNAR)|54.1|46.2|32.9|23.9|\n|Mixed (MCAR+MAR+MNAR)|34.9|23.8|25.5|22.8|\n\nIn summary, the dependency of our algorithm's complexity on the degree of missingness is an inverse relationship. Practitioners can expect our method to run faster on datasets with higher proportions of missing values, which we assume is a desirable property for real-world applications.\n\n## **References**\n[1] Kwon, Jeongyeol, and Constantine Caramanis. \"The EM algorithm gives sample-optimality for learning mixtures of well-separated gaussians.\" Conference on Learning Theory. PMLR, 2020.\n\n[2] Celeux, G., Maugis-Rabusseau, C., & Sedki, M. (2019). Variable selection in model-based clustering and discriminant analysis with a regularization approach. Advances in Data Analysis and Classification, 13, 259-278. Springer.\n\n[3] Maugis, C., Celeux, G., & Martin-Magniette, M-L. (2009). Variable selection in model-based clustering: A general variable role modeling. Computational Statistics & Data Analysis, 53(11), 3872-3882. Elsevier.\n\n[4] Witten, D. M., Friedman, J. H., & Simon, N. (2011). New Insights and Faster Computations for the Graphical Lasso. Journal of Computational and Graphical Statistics, 20(4), 892-900.\n\n[5] Mazumder, R., & Hastie, T. (2012). Exact Covariance Thresholding into Connected Components for Large-Scale Graphical Lasso. Journal of Machine Learning Research, 13(27), 781-794."}, {"review_id": "mQ1Vq6QIuc", "reviewer": "Reviewer_Wfhz", "summary": "This paper introduces a unified framework for variable selection in model-based clustering when data contain Missing Not at Random (MNAR) entries. The key contribution is the development of a penalized likelihood approach using a joint model for both the data and the missingn mechanism. The authors propose an efficient EM algorithm for estimation and selection and demonstrate the method’s effectiveness through simulation studies and real data applications.", "questions": "Please see the Weaknesses.", "limitations": "Please see the Weaknesses.", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n\nThe paper focuses on variable selection in model-based clustering under MNAR, which is rarely addressed. The proposed method is computationally efficient and theoretically well-justified, relying on an EM-type iterative procedure specifically adapted to handle MNAR data structures. The experimental design is thoughtful and comprehensive, effectively demonstrating the superior performance of the proposed method across various missing scenarios. Furthermore, the consistency of the variable selection procedure is rigorously established under suitable regularity conditions.\n\nWeaknesses:\n\n1. One of my primary concerns regarding this manuscript is the level of novelty, as the proposed method appears to be an extension of [1] and [2]. The authors need to justify the novelty of the proposed method compared to these works, from either methodological or theoretical perspective. Therefore, I hold a conservative opinion toward the significance of the contribution made by this manuscript.\n2. The paper assumes that the missing mechanism follows parametric regression model. This can be restrictive in complex practical cases.\n3. There is no numerical studies or test of model misspecification for missing mechanism models.\n4. The method is currently limited to continuous data modelled by Gaussian mixtures. Many practical applications (e.g., survey or electronic health record data) involve categorical or mixed-type variables. This limitation is briefly acknowledged but not substantively addressed.\n5. While the paper claims that the proposed EM-based framework is computationally feasible, this claim is not supported by quantitative analysis or complexity evaluation. This is a notable limitation, particularly given the well-documented computational challenges associated with EM algorithms in high-dimensional settings or when the number of clusters is large.\n6. This paperdoes not report runtime, memory usage, or scalability evaluations, making it difficult to assess the method’s practicality on large-scale datasets. Including a detailed assessment of computational performance, such as runtime comparisons across competing methods, scalability experiments with increasing dimensionality and sample size, or at least a discussion of algorithmic complexity, would substantially enhance the empirical contribution.\n7. The method depends on several hyper-parameters (e.g., $\\lambda$, $\\rho$, and the stopping rule parameter $c$), yet no sensitivity analysis or guidance is provided for their selection.\n\n[1] Maugis-Rabusseau, C., Martin-Magniette, M. L., & Pelletier, S. (2012). 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.\n\n[2] Celeux, G., Maugis-Rabusseau, C., & Sedki, M. (2019). Variable selection in model-based clustering and discriminant analysis with a regularization approach. Advances in Data Analysis and Classification, 13, 259-278.", "was_revised": true, "final_justification": "While I appreciate the authors’ responses, I maintain my original evaluation, as the contribution remains relatively limited and the literature review is insufficient.", "rebuttal": "We thank the reviewer for their time and for providing a detailed and structured critique of our manuscript. The feedback is very valuable, and we appreciate the opportunity to clarify our contributions and address the perceived limitations. While the reviewer raises several important points, we believe some of the primary concerns, particularly regarding novelty and the handling of MNAR data, stem from a key distinction between our work and the cited references. We will address each weakness point-by-point.\n\n### **1. On Novelty Compared to Maugis et al. [1] and Celeux et al. [2]**\n\nWe respectfully but strongly disagree with the reviewer's assessment. Your concern is understandable, as we build upon the SRUW framework and SelvarMix that Maugis and Celeux pioneered. However, the cited papers operate under a fundamentally different and simpler assumption about the missing data, which lies in the critical distinction between MAR and MNAR:\n\n+ Maugis et al. [1] (SelvarClustMV) explicitly and exclusively handle data that is **Missing at Random (MAR)**. Their EM algorithm is designed for the MAR case, where the missingness mechanism can be safely ignored after conditioning on observed data.\n\n+ Celeux et al. [2] (SelvarMix) propose a two-stage penalized approach that also operates within the standard likelihood framework, which is only valid for MAR or MCAR data.\n\nOur work is the first to extend this entire variable selection paradigm to the much more challenging Missing Not at Random (MNAR) setting. The reasons this is not a trivial extension are:\n\n1.  **Methodological aspect:** Moving from MAR to MNAR is not a simple extension. It requires developing a **new joint model** for the data and the missingness mechanism (our SRUW-MNARz model) and modifying a **new EM algorithm** to handle the non-ignorable nature of the missingness. As we have presented in the Appendix, the GMM representation of SRUW under MAR, its identifiability, and the equivalence of MAR and SRUW under MNARz are essential for theoretical guarantee and modeling; thereby, they are novel. The likelihood function, the E-step, and the M-step are all fundamentally different from those in [1] and [2].\n\n2.  **Theoretical aspect:** The theoretical guarantees for identifiability and selection consistency provided in our **Theorems 1, 2, and 3** are entirely new and, in some sense, complement Celeux's work as it concurrently serves as one of our objectives. The proofs in the cited works are not directly applicable to the MNAR setting, as the statistical properties of the joint likelihood are far more complex. Establishing these guarantees for our new model class is also a core theoretical contribution of our paper.\n\n### **2. On the Parametric Missingness Mechanism and Model Misspecification**\nWe agree that you have a valid point. Our choice of the MNARz mechanism is a deliberate trade-off between **model flexibility** and **statistical identifiability**.\n\n+ As shown by recent literature (e.g., Sportisse et al. [3]), general non-parametric MNAR models are often non-identifiable without external information, making inference impossible. The parametric MNARz model is one of the few known mechanisms that is **provably identifiable** in a clustering context, allowing for consistent parameter estimation. We believe that starting with a theoretically sound, identifiable model is a more rigorous scientific approach than using a more flexible but potentially unidentifiable one.\n\n+ We acknowledge that we did not include a formal misspecification test. However, we did perform an empirical robustness check. In **Appendix E.2 (Figure 4)**, we tested our method on data generated under `MNARy` and `MNARyz` mechanisms similar to Sportisse et al. [3], where the missingness depends on the unobserved values themselves, a direct violation of our model's assumption. Furthermore, as suggested by reviewer BwuH, we also conducted a study under the case mixture of missing mechanisms. Our method still performed competitively and often outperformed baselines, demonstrating a degree of robustness to this form of misspecification.\n\n### **3. On Limitations (Continuous Data Only)**\nWe fully agree with the reviewer on this point. This is a current limitation of our framework. Extending model-based clustering with variable selection and MNAR handling to mixed-type data is a highly complex research program in its own right. Our work provides the foundational framework for continuous data, which is a necessary and significant first step. However, we believe this framework can be extended naturally, and we assure you that this is not just a hypothetical future direction but an active area of research. For instance:\n+ An extension to general mixed data would likely follow a similar path, integrating ideas from model-based approaches, such as Marbac & Sedki [4]. It is crucial to note, however, that their framework is designed for the simpler MAR case. Adapting such a model to handle the non-ignorable MNAR mechanisms, as we do here, would be a significant research effort in its own right. Moreover, the variable roles in their model, where all variables are pairwise independent, are also much simpler than the current SRUW that we adopted in our framework.\n\n+ Jacques & Murphy [5] successfully applied this exact philosophy to multivariate count data; however, in a simpler setting where they only have two variable roles rather than three, like the one we use. The fundamental logic remains the same, but they replaced the models for continuous data with Poisson mixtures and Poisson Generalized Linear Models to correctly handle the discrete nature of the variables. This shows that the underlying variable selection strategy is not tied to a single data type.\n\nWe will revise the discussion to include this context, clarifying both the current scope and these evidence-based pathways for future extensions, placing it within a broader research trajectory.\n\n### **4. On Lack of Computational Analysis and Hyperparameter Guidance**\nWe thank the reviewer for this crucial feedback and agree that this was a significant omission in our original submission. To rectify this, we have conducted an informal computational analysis, which we will add to the supplementary material.\n\nOn the model's complexity, we respectfully agree that while the model is complex, the problem we address in this paper is inherently complicated. Simple models, like the impute-then-cluster baselines, fail precisely because they do not capture the joint structure of the data, variables, and missingness mechanism.\nOur framework's complexity is a necessary consequence of its principled and comprehensive approach.\n\nOn hyperparameter tuning, we agree that more guidance would be beneficial.\n\n1. For the stopping criterion, this parameter balances false positives and false negatives in the variable selection step and we also provided a sensitivity analysis in Appendix (Appendix E.2). Our results support the previous finding by Celeux et al. [1] that the value of c between 3 and 5 provides a good trade-off, and performance is stable within this range for well-posed problems\n\n2. For the regularization parameters $\\lambda$ and $\\rho$, they can be selected via a grid search, which is a standard approach for penalized methods. In our experiments, we do not perform an arbitrary search but rather a principled and data-driven procedure to automatically construct a useful grid, which ensures the search space is relevant to the data at hand. The process involves these main steps:\n  + Derive $\\lambda_{\\text{max}}$ from the maximum initial score for the mean parameters, $\\text{max}|t(Z_0) X|$.\n  + Derive $\\rho_{\\text{max}}$ as a similar upper bound for the precision matrix penalty, derived from the initial sample covariances and cluster responsibilities.\n  + The two values above are motivated to be the smallest penalties that would shrink all corresponding parameters to zero.\n  + We then generate a geometric sequence of $L$ values (e.g., $L = 5$) from a small fraction of the maximum value up to the maximum value itself. This standard practice ensures that the grid efficiently explores penalties across several orders of magnitude.\n\nRegarding the scalability problem, we appreciate the reviewers' insightful feedback and their important questions regarding our framework's computational performance. The concerns about complexity, scalability, and practical feasibility in high-dimensional settings are indeed valid and critical for this method. To address this, we've conducted a detailed study, combining a sketchy, informal treatment of complexity with an experimental analysis. Our results demonstrate that our framework is not only theoretically favorable but also computationally practical, offering a polynomially scaling speedup, which makes it a capable tool for the problems it's designed to solve. The more detailed analysis can be found in the rebuttal of reviewer rfVw. \n\n### **References**\n[1] Celeux, G., Maugis-Rabusseau, C., & Sedki, M. (2019). Variable selection in model-based clustering and discriminant analysis with a regularization approach. Advances in Data Analysis and Classification, 13, 259-278. Springer.\n\n[2] Maugis, C., Celeux, G., & Martin-Magniette, M-L. (2009). Variable selection in model-based clustering: A general variable role modeling. Computational Statistics & Data Analysis, 53(11), 3872-3882. Elsevier.\n\n[3] A. Sportisse, M. Marbac, F. Laporte, G. Celeux, C. Boyer, J. Josse, and C. Biernacki. Model\nbased clustering with missing not at random data. Statistics and Computing, 34(4):135, 2024.\n\n[4] M.Marbac and M.Sedki. Variable selection for mixed data clustering: A model-based approach. arXiv preprint arXiv:1703.02293, 2017.\n\n[5] Jacques, J., & Murphy, T. B. (2025). Model-Based Clustering and Variable Selection for Multivariate Count Data. Computo."}, {"review_id": "MnchWid9F1", "reviewer": "Reviewer_Bx5V", "summary": "This paper introduces a model for variable selection in finite Gaussian mixtures. The basic idea is to partition the data vector into variables that are relevant for the clustering, variables that are linear combinations of the variables relevant to the clustering, and then irrelevant variables. To accomodate missing not at random variables, the authors append a missingness mechanism to the finite mixture model, which they call MNARz, that allows for missingness probabilities to vary within a latent class. They then use a penalized likelihood function to perform simultaneous parameter estimation and variable selection.\n\nThe authors apply their model for clustering, variable selection, and parameter inference in a simulation study, where the method demonstrates compentitive imputation accuracy and clustering performance. In their real data application, they use their model to analyze transcriptome data, where they determine scientifically relevant subgroups in the data based on gene expressions.", "questions": "1) Estimation and selection: I am quite confused about how variables are simultaneously selected and the finite GMM is estimated. The paper would benefit from clearer and higher-level algorithms describing estimation and selection. From the Finite GMM representation of the SRUW model, does a user need first to declare which variables they believe are relevant/irrelevant to the clustering? When/Where is the penalty applied? There is some mention of a ranking step followed by a parameter estimation step but this is not clear in the paper or manuscript. Where/when do the penalties come into play?\n\n2) MNARz mechanism: I have some serious doubts that the proposed MNARz model provides a missingness mechnism that is nonignorable. This needs to be clarified. For example the authors state that the missignness mechanism within any class factors as:\n\np(c_{n} \\mid z_{n_k}, \\psi_{k}) = \\prod_{d=1}^{D}\\rho(\\psi_{k})^{c_{nd}}(1-\\rho(\\psi_{k})^{1-c_{n})\n\nFor the missingness to be MNAR, missingness has to depend on unobserved missing values. In this case, this would require \\psi_{k} to depend on y_{d}, but this is not clear from how the model is presented. Also, the way this reads is that \\psi_{k} is the same for each variable y_d in latent class k. Is this correct? This seems highly improbably and is not a MNAR missigness mechanism. This needs to be clarified or modified if the authors claim their methods are useful for MNAR data. \n\n\n3) The distinction between the different \\pi_{k}s in (3) is not clear.", "limitations": "1) In general, the paper is heavy with notation that is difficult to follow\n2) Following my above questions:\n- I think there should be some outline for how model estimation and variable selection is conducted. The E-M algorithm presented in the main text only focuses on estimating the parameters in the joint model for the missingness and study variables for a fixed set of pertinent/irrelevant variables, unless I am missing something\n-I think the authors should revisit the MNARz missingness mechanism. In general, a variable is nonignorably missing if the probability of missingness depends on the underlying missing values. I don't see how this is encoded in the MNARz variables\n-Related, I also searched in the supplement for how the missingness mask is realized but couldn't find anything in the text or supplement. Is the data truly MNAR in the author's simulated experiments?", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 1, "significance": 3, "originality": 3, "strengths_and_weaknesses": "I think the strength of this paper is in the ability of the SRUW model to be expressed as a straightforward finite Gaussian mixture, with component-wise mean vector and covariance matrix expressed in terms of the relevancy partition that the authors originally present. As the authors mention, this enables theoretical analyses and relatively simple E-M algorithms for estimation. In addition, the idea that clustering is driven by a few of the study variables, and seeking to select these variables to provide interpretable insights, is quite attractive. Extending these ideas to MNAR study variables is quite challenging.", "was_revised": true, "final_justification": "Once the authors revise their work according to their detailed responses, the paper will be satisfactory. For this reason, I recommend boarderline accept.", "rebuttal": "We sincerely thank you for your thoughtful and detailed feedback. The questions raised are extremely helpful and have highlighted several areas where our presentation was not as clear as it should have been. We agree that the paper is dense with notation and some parts of the writing are unclear, and we are grateful for the opportunity to clarify our methodology, particularly the variable selection procedure and the nature of the MNARz mechanism. We will address each of your points below.\n### **1. On Clarity of Estimation and Variable Selection**\nWe deeply apologize for the lack of clarity here. Your confusion is our responsibility, and we will significantly revise Section 3 and add a high-level algorithm to make this process transparent. To answer the questions directly:\n+ No, the user does not pre-declare variable roles. The entire purpose of our framework is to discover these roles (relevant, redundant, irrelevant) from the data automatically. However, if the user suspects that the data can include both MNAR and MAR mechanisms, they have to provide the variable roles beforehand to use the mixed mechanism modelling, which is currently our limitation on this mixed mechanisms handling part.\n+ The penalties are applied in the first stage (Variable Ranking). This is the key to making the problem computationally tractable since the result will be a list of ranked variables and we just have to do a **linear pass**\n+ The process is a two-stage procedure, as pioneered by Celeux et al. [1], which we have extended to the missing data setting.\nHere is the high-level description of the algorithm, which we will add to the revised manuscript:\n\n**Stage 1: Variable Ranking via Penalized GMM**\n+ **Target:** To obtain an ordered list of all $D$ variables, from most to least relevant for clustering.\n+ **Method:** We fit a penalized Gaussian Mixture Model (GMM) to the complete data $X$ where missing entries in $X$ has been imputed naively. The reason to perform this step is that we want to have a good ranking scores and it is critical to note that this imputed dataset is only a temporary scaffold used to enable the ranking procedure; it is discarded immediately after.\nThe ranking itself is then performed by fitting a penalized log-likelihood given by Equation (3) in the manuscript to the temporary complete data. The $\\lambda$ penalty shrinks the cluster means $\\mu_k$ towards a common center. Variables that are important for clustering will resist this shrinkage. The $\\rho$ penalty, weighted by the adaptive matrices $P_k$, shrinks the off-diagonal elements of the precision matrices $\\Psi_k$ to zero, effectively performing graphical model selection within clusters.\n\n+ **Procedure:** We fit this model over a grid of $(\\lambda, \\rho)$ values. For each variable d, we compute a \"ranking score\" $\\mathcal{O}\\_{K}(d)$. This score quantifies how essential the variable is for defining the cluster structure across the range of penalties. In essence, a variable that is important for clustering will exhibit significant separation in its cluster-specific mean estimates $\\hat{\\mu}\\_{kd}$ even when the regularization penalty is strong. The score $\\mathcal{O}\\_{K}(d)$ is a summary statistic that captures this behavior.\n+ **Output:** A single ranked list of variables $(d_1, d_2, ..., d_D)$.\n\n**Stage 2: Role Assignment using the SRUW Model and BIC**\n+ **Target:** To partition the ranked list into the sets $\\mathbb{S}$ (relevant), $\\mathbb{U}$ (redundant), and $\\mathbb{W}$ (irrelevant).\n+ **Method:** We perform a single, efficient linear pass through the ranked list of variables.\nWe start with a small model and iteratively add the next variable from the list.\nAt each step, we fit the full (unpenalized) SRUW-MNARz model and use the BIC criterion (Equation (2)) to decide the role of the newly added variable (i.e., whether it belongs to  $\\mathbb{S}$, $\\mathbb{U}$, and $\\mathbb{W}$).\nThis avoids the computationally prohibitive $\\mathcal{O}(D^5)$ complexity of the original stepwise SRUW method by replacing its combinatorial search with a single pass over a statistically-informed ordering. Crucially, this two-stage procedure is precisely what enables the \"relatively simple E-M algorithms\" and theoretical construction that you correctly identified as a key strength of our work.\n\n### **2. On the MNARz Mechanism and its Non-Ignorable Nature**\nYou are correct in your understanding of the classical definition of MNAR. However, the MNARz mechanism is a specific, identifiable subclass of MNAR models that has been recently established in the statistical literature.\n\n#### **2.1. Why MNARz is Truly MNAR (Non-Ignorable):**\nIn a mixture model context, the latent variable $z_n$ (the cluster assignment) is itself part of the unobserved data. The missingness probability $P(c_{nd}=1 | y_n, z_n)$ depends on $z_n$. Since $z_n$ is unobserved and must be integrated out, the marginal probability of missingness $P(c_{nd}=1 | y_n^{obs})$ will depend on the posterior distribution $P(z_n | y_n^{obs})$, which in turn depends on the observed data $y_n^{obs}$. Crucially, because the parameters of the data distribution $(\\theta_k)$ and the missingness distribution $(\\psi_k)$ are linked through the same latent variable $z_n$, they cannot be estimated separately. **This linkage** is what makes the mechanism **non-ignorable**. If one were to ignore the missingness mechanism, the estimation of the mixture parameters $\\theta_k$ and the cluster assignments would be biased. This is the formal definition of a non-ignorable mechanism.\nThis model class was formally introduced and its properties studied by Sportisse et al. (2024) (see [1]). It represents a principled and tractable way to move beyond MAR when the cause of missingness is related to the underlying latent structure (e.g., in transcriptomics, genes in a specific functional cluster might fail to express and thus be \"missing\" due to that cluster's biological state).\n\n#### **2.2. Why $(\\psi_k)$ is Constant Across Variables Within a Class**\nYou are correct to point out that in our current formulation, the missingness probability $\\rho(\\psi_k)$ is a single scalar for all variables $d$ within a given class $k$. This is a **deliberate** modeling choice made for two reasons:\n\n*   **Parsimony and Estimation Stability:** In a high-dimensional setting ($D \\gg N$), introducing a separate missingness parameter $\\rho_{kd}$ for each variable in each cluster would add $KD$ parameters. This would make the missingness mechanism highly over-parameterized and prone to severe overfitting, leading to unstable estimation. Our approach, by contrast, adds only $K$ parameters. This deliberate choice for parsimony ensures the model remains statistically identifiable and that the EM algorithm converges to a stable solution, preventing the core clustering structure from being obscured by noise.\n\n*   **Enhanced Interpretability:** This choice leads to a highly interpretable result. A high value of $\\rho_k$ for a specific cluster $k$ provides a simple, powerful insight: \"Cluster $k$ represents a subgroup (e.g., a biological state or patient type) that is systematically harder to measure or has lower quality data across the board.\" This single parameter per cluster is far more interpretable than a dense matrix of $KD$ individual probabilities. Ultimately, this modeling choice for parsimony is how we make the problem of MNAR modeling stable and tractable, directly addressing the challenge you highlighted while preserving the interpretability you found attractive.\n\nMoreover, in our simulations, the missingness mask is generated precisely according to this class-dependent MNARz mechanism, ensuring the simulated data truly reflects this non-ignorable pattern.\n\n### **3. On the Distinction Between $\\pi_k$ in Equation (3)**\n\nWe apologize for this confusing notation, and this is a clear notational flaw on our part when creating the manuscript. They are actually $\\pi_k$ and $P_k$. Here is the clear explanation:\n+ $\\pi_k$ refers to the mixing proportions of the GMM, i.e., $P(z_n=k)$. It is the scalar mixing proportion for component $k$.\n+ $P_k$ in the penalty term $(\\mathbf P_K \\odot \\mathbf\\Sigma_k^{-1})_{dd'}$ refers to the adaptive weight matrix used for the graphical Lasso penalty. As we mentioned briefly, these are typically chosen to be inversely proportional to initial partial correlations, following the adaptive Lasso principle.\n\n### **References**\n\n[1] A. Sportisse, M. Marbac, F. Laporte, G. Celeux, C. Boyer, J. Josse, and C. Biernacki. Model\nbased clustering with missing not at random data. Statistics and Computing, 34(4):135, 2024."}]

Author rebuttal:

Meta-review: The authors introduce a method for selecting variables in model-based clustering that effectively addresses missing values, even when they are not randomly distributed. The solution is based on a data-driven penalty and explicitly models the missing data through latent class membership. The method uses an extended EM algorithm, and the authors provide theoretical guarantees of consistency of the solution. Additionally, it shows empirical enhancements in the identification of informative variables and the clustering of complex datasets. The problem is particularly relevant in bioinformatics, where missing values are prevalent and variables are often not informative. Most reviewers were positive about this paper and acknowledged the methodological and theoretical contribution of this work. The discussion was active, and most issues were addressed by the authors. Specifically, the authors provided additional results supporting the effectiveness of the method. They also implemented another two-stage approach as a baseline; the results support the advantage of the proposed optimization problem. One reviewer had a remaining concern regarding novelty. Still, after looking at the related work, the work differs enough from existing results that it should be considered a valuable non-trivial contribution. 

Overall, I support accepting this paper. While the problem is a bit niche, this is a valuable new idea that leads to empirical benefits, stability, and is theoretically justified. Additionally, the paper is well written, and the approach could find applications in bioinformatics.
Decision: Accept (poster)