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

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.

1. On Clarity of Estimation and Variable Selection

We 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:

• 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.
• 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
• The process is a two-stage procedure, as pioneered by Celeux et al. [1], which we have extended to the missing data setting. Here is the high-level description of the algorithm, which we will add to the revised manuscript:

Stage 1: Variable Ranking via Penalized GMM

• Target: To obtain an ordered list of all $D$ variables, from most to least relevant for clustering.

• 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. The 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.

• 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.

• Output: A single ranked list of variables $(d_1, d_2, ..., d_D)$.

Stage 2: Role Assignment using the SRUW Model and BIC

• Target: To partition the ranked list into the sets $\mathbb{S}$ (relevant), $\mathbb{U}$ (redundant), and $\mathbb{W}$ (irrelevant).
• Method: We perform a single, efficient linear pass through the ranked list of variables. We start with a small model and iteratively add the next variable from the list. At 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}$). This 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.

2. On the MNARz Mechanism and its Non-Ignorable Nature

You 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.

2.1. Why MNARz is Truly MNAR (Non-Ignorable):

In 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. This 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).

2.2. Why $(\psi_k)$ is Constant Across Variables Within a Class

You 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:

• 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.

• 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.

Moreover, 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.

3. On the Distinction Between $\pi_k$ in Equation (3)

We 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:

• $\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$.
• $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.

References

[1] A. Sportisse, M. Marbac, F. Laporte, G. Celeux, C. Boyer, J. Josse, and C. Biernacki. Model based clustering with missing not at random data. Statistics and Computing, 34(4):135, 2024.

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.

1. On Novelty Compared to Maugis et al. [1] and Celeux et al. [2]

We 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:

• 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.

• 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.

Our 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:

1. 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].

2. 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.

2. On the Parametric Missingness Mechanism and Model Misspecification

We agree that you have a valid point. Our choice of the MNARz mechanism is a deliberate trade-off between model flexibility and statistical identifiability.

• 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.

• 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.

3. On Limitations (Continuous Data Only)

We 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:

• 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.

• 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.

We 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.

4. On Lack of Computational Analysis and Hyperparameter Guidance

We 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.

On 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. Our framework's complexity is a necessary consequence of its principled and comprehensive approach.

On hyperparameter tuning, we agree that more guidance would be beneficial.

1. 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

2. 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:

• Derive $\lambda_{\text{max}}$ from the maximum initial score for the mean parameters, $\text{max}|t(Z_0) X|$.
• Derive $\rho_{\text{max}}$ as a similar upper bound for the precision matrix penalty, derived from the initial sample covariances and cluster responsibilities.
• The two values above are motivated to be the smallest penalties that would shrink all corresponding parameters to zero.
• 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.

Regarding 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.

References

[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.

[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.

[3] A. Sportisse, M. Marbac, F. Laporte, G. Celeux, C. Boyer, J. Josse, and C. Biernacki. Model based clustering with missing not at random data. Statistics and Computing, 34(4):135, 2024.

[4] M.Marbac and M.Sedki. Variable selection for mixed data clustering: A model-based approach. arXiv preprint arXiv:1703.02293, 2017.

[5] Jacques, J., & Murphy, T. B. (2025). Model-Based Clustering and Variable Selection for Multivariate Count Data. Computo.

Thank you for the detailed follow-up and for providing these specific references. We will clarify the scope of our work and its novelty.

1. On the Novelty of Variable Selection under MNAR

We appreciate you raising these important papers. We agree that variable selection under MNAR is an active research area. However, we believe there is a fundamental distinction between the problem these papers solve and the one we address.

All the cited works focus on supervised learning (regression) settings, where the goal is to select important predictors for a known outcome variable. Even in complex cases like Mixture of Experts regression, the clustering is guided by how well different models predict this known outcome.

Our work, by contrast, addresses the unsupervised learning task of model-based clustering. Here, there is no outcome variable; the goal is to simultaneously discover the unknown group structure (the clusters) and identify the variables that define that structure. This is a fundamentally different problem, as the variable roles (relevant, redundant, irrelevant) are defined with respect to a latent structure that must be learned from the data itself.

Therefore, we will refine our claim to be more precise: "To our knowledge, this is the first work to extend the SRUW variable-role modeling framework to the MNAR setting for model-based clustering." This distinguishes our contribution from the important work in the supervised context.

2. On Model Misspecification

Regarding the issue of model misspecification, we apologize if this was not clear in our initial rebuttal. We have indeed conducted these studies, and the results are in the submitted supplementary material.

• Existing Robustness Check: Appendix E.2 (Figure 4) contains a robustness check where we test our method on data generated under MNARy and MNARyz mechanisms, which are direct violations of our model's assumption.

  • MNARy: where missingness depends directly on the unobserved value itself, independent of the cluster.

  • MNARyz: a more complex case where missingness depends on both the unobserved value and the latent cluster membership.

• New Mixed-Mechanism Experiments: Furthermore, in response to the reviewer BwuH, we conducted new experiments on mixed-mechanism data (MCAR+MAR+MNAR), which we will add to the appendix.

In all these studies, our method consistently outperforms MAR-based approaches, demonstrating its practical robustness.

3. On the Sensitivity Analysis

You are absolutely right to have been unable to find the results, and we sincerely apologize for our error. In our previous response, we mistyped the location. The sensitivity analysis is located in Appendix E.3, not E.2. While we have already provided the sensitivity results for the stopping parameter c, you are correct that we did not run a new simulation for the regularization parameters $\lambda, \rho$.

We would like to clarify that this was a deliberate choice based on the observed properties of the two-stage, regularization-path-based ranking method, which was pioneered by Celeux et al. (2019). The robustness of the variable ranking to the specific choice of the $\lambda, \rho$ grid is a key feature of this approach.

The ranking score, $\mathcal{O}_{K}(d)$, is not based on a single model fit at an arbitrary $(\lambda, \rho)$ value. Instead, it summarizes the behavior of a variable across the entire regularization path. Irrelevant variables have their parameters shrunk to zero very early and stay there, while relevant variables resist this shrinkage. This makes the final ranking quite stable to the precise granularity of the grid.

This exact point was discussed by the originators of the method. Here, we summarize their findings:

1. They state that the ranking criterion "is weakly dependent on the chosen grid" because irrelevant variables are detected early in the path.
2. They also note that for this ranking stage, they deliberately ignore the variance matrices to reduce complexity, justifying it because the primary interest is in separating clusters by their means.

Our work follows their justified methodology. We use a principled, data-driven method to construct the grid (as described in our rebuttal), and we also pre-center and scale our data to align with their second remark.

Therefore, while a dedicated simulation for $(\lambda, \rho)$ is possible, the known robustness of the ranking score itself led us to believe that focusing our new empirical studies on the more critical and previously unanswered question of misspecification of the missingness mechanism was a more significant contribution to validating our framework.

We will add a paragraph to the appendix explicitly discussing this point and referencing the prior work to provide this important context for the reader. We hope this clarifies our reasoning and addresses your concern. Thank you again for your constructive engagement.

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.

1. On the Dependence on the MNARz Scheme and Robustness

Theoretically, 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.

In 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:

1. MNARy: The missingness probability depends on the unobserved value $y_{nd}$ itself.

2. MNARyz: The missingness probability depends on both the unobserved value $y_{nd}$ and the latent class $z_n$.

We 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.

2. On the Initialization of the EM Algorithm

EM for finite mixtures is famously sensitive to initialization. So, we use a standard and effective multi-step initialization procedure:

1. 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.
2. We apply hierarchical agglomerative clustering or $K$-means on this completed dataset to obtain an initial partition of the data into $K$ clusters.
3. 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.

This 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.

Table 4: Comparison of EM Initialization Methods under 8 Scenarios

MIs: Multiple Initializations

3. On Computational Complexity and Scalability

We thank the reviewer for their crucial feedback. We 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:

1. 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$.

2. 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.

3. 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.

4. 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.

3.1 Classical Backward Stepwise SRUW

This 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.

• Fitting a model with $d$ variables costs $\mathcal{O}(M_{\text{EM}}\cdot (NKd^2+Kd^3))$.

• 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.

• The total removal cost is:

$$C_{\text{stepwise}}=\mathcal{O}(M_{\text{EM}}\cdot (NKD^4+KD^5))$$

The $\mathcal{O}(D^5)$ complexity makes this approach computationally infeasible for even moderately high dimensions (e.g., $D>100$).

3.2 Two-Stage Lasso-Like SRUW

Our proposed method avoids this combinatorial explosion.

• 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]).

$$C_{\text{rank}}=\mathcal{O}(M_{\text{grid}}\cdot M_{\text{EM}}\cdot(NKD^2+KM_{\text{glasso}}D^2))$$

• 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$.

Under the sparsity assumption ($D_{\text{eff}}\ll D$), the overall complexity is: $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}}))$

3.3 Computational Speedup

The speedup is the ratio $\frac{C_{\text{stepwise}}}{C_{\text{lasso}}}$.

• 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}}})$.
• 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}}})$.

In 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.

Table 5: Empirical Runtime (in seconds) and Speedup Factor

4. On the Dependency of Complexity on the Degree of Missingness

We have conducted a new empirical study to provide a definitive answer.

Our 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:

Table 6: Empirical Runtime (in seconds) with Increasing Missing Rate under Mixed Missingness Mechanism

In 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.

References

[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.

[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.

[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.

[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.

[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.

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

1. On handling a mixture of missing data types

Our 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:

• 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.
• 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.
• We denote "MNARz + SelvarMix" to be the suggested baseline from the reviewer. More details are in answer for question 3.
• 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)
• 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.

Here is the result:

Table 1: Performance on Mixed (MAR+MNARy) Data

Table 2: Performance on Mixed (MCAR + MAR + MNARy) Data

2. On the Clarity and Interpretation of Figure 1

We 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:

• 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.
• 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.
• 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.

Table 3: Welch's t-test for Mean ARI between SelvarMNARz and Clustvarsel with $\alpha = 0.05$

• 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.

3. On Decoupled versus Joint Approaches

We 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":

• 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.
• 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.

The 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.

• 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.
• 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.