$\text{G}^2\text{M}$: A Generalized Gaussian Mirror Method to Boost Feature Selection Power

Thank you for your time and feedback, and for assessing that the theoretical claims are easy to follow and numerical experiments are extensive and detailed.

QB1: With a rich background of the feature selection problem and extensive existing methods mentioned by the authors, the paper does not formulate or describe any existing method in the main text. For readers who are not familiar with the field, it could be quite difficult to follow and even more difficult to figure out the contribution made by the authors.

We put the two related works– the Gaussian mirror and data splitting methods in Appendix A. And this was primarily due to the page limit. To improve clarity, we will move Appendix A to the main paper in the revised manuscript..

QB2: In comparison, the details for simulation settings and implementation take up too much space. I believe it is the main reason for weakness 1 since the page limit is 9. I strongly recommend that the authors delay unnecessary details of the experiments to the appendix and save more space for the introduction part.

Thank you for the suggestion. We will shorten the description of each of the datasets used in the experiment section in the revised manuscript and also move the full description to a separate section in the Appendix to balance the introduction of the related work section according to QB1.

QB3: Since the authors assume a linear model between X and Y, is G$^2$M valid for general regression models? Or is it possible to extend G$^2$M to general models?

It is possible to extend G$^2$M to the generalized linear model setting. For example, Reference [1] discusses the generalized linear model setting with Gaussian mirror statistics. We believe there is a similar extension for G$^2$M. However, given that the form of test statistics in our paper is significantly different from the original Gaussian mirror, the extension of G$^2$M to this case is non-trivial. And this is deferred to future work. We will add the above discussion in Section 4 of the revised manuscript.

QB4: The coefficients $\beta^+$ and $\beta^-$ are used without clear definitions. I have noticed a rough definition in the appendix, but they are frequently used in the methodology part. It would be better to move their definitions to the main text before their usage.

According to the response to QB1, we moved the related work section from Appendix A to the main paper in the revised manuscript. This should provide enough context on how the coefficients $\beta^+$ and $\beta^-$ are generated from the data splitting and the Gaussian mirror context. Essentially, the $\beta^+$ and $\beta^-$ used in G$^2$M follow the same generation process in a variance generalized setting. We will explicitly mention this connection in the updated related work section for the Gaussian mirror in the revised manuscript.

QB5: Although the method is based on existing works that ensure FDR control, it would be better to explicitly formulate this property in a short theorem. It might be a straightforward result or even only a proposition, but showing this issue is worthwhile since there are also many intuitive methods for feature selection that have no error rate control guarantee.

Thanks for the suggestion. In the revised manuscript, we will state the theorem on FDR control for data splitting in the revised manuscript. Specifically. We state:

Theorem 1: For a given set of mirror statistics $w_j$, $j \in [p]$, if mirror statistics satisfy the following property: 1. If $j \in S_0$ is a null feature index, then given any real number $t \ge 0$, $p(w_j \ge t) = p(w_j \le -t) = 0.5$. 2. If $j \in S_1$ is a non-null feature index, then $p(sign(w_j) = 1) = 1$.

Then given a nominal FDR level $q$, the feature selection set $\lbrace j \vert w_j \ge \tau_q\rbrace$ controls FDR, where $\tau_q$ is chosen via Eq.(3).

The proof idea is already informally presented in lines 469-473, and we will provide the following proof for this Theorem in the revised manuscript.

Proof: Given the above two properties for the mirror statistics, we realize that $\vert \lbrace j: w_j \le -t, j\in S_0\rbrace \vert = \vert \lbrace j: w_j \ge -t, j\in S_0\rbrace \vert$ due to the symmetric-about-zero property for the mirror statistics $w_j$ under the null distribution. As a result, the selection rule $\hat{S}_1 = \\lbrace w_j \ge \tau_q\rbrace$ (e.g., the denominator of Eq. (3)) selects features based on the level-$q$ quantile $\tau_q$ on the right side of the null distribution, ensuring the FDR level to be at most $q$.

References

[1] Dai, C., Lin, B., Xing, X., & Liu, J. S. (2023). A scale-free approach for false discovery rate control in generalized linear models. Journal of the American Statistical Association, 118(543), 1551-1565.

Thank you for your time and feedback, and for assessing that the paper has a strong theoretical foundation and the methodology provides greater empirical power compared to benchmarks.

QA1: The paper does not discuss whether the proposed test statistics have an intuitive interpretation, which may be of interest in the case of real-world data.

Essentially, the proposed test statistics is mainly the result of a variance extension of the mirror statistics. We have included in Appendix B why it is necessary to consider a variance-generalized setting. Primarily, this is because we observed that even with the simplest Gaussian simulation setting (see Appendix B), the variance of the mirror statistics is not 1. When data is more complex (e.g., the real-world data), it is expected that the unit-variance setting will be violated. In this case, the original Gaussian Mirror/ Data splitting setting that depends on the unit-variance assumption is likely to fail. In QA3 below, we have also included another real-world dataset and demonstrate the gains of the proposed statistics from the feature selection perspective.

QA2: There is still some ambiguity concerning the difference between G$^2$M and G$^2$M$^\dagger$. Given that the highest statistical power is predominantly achieved using G$^2$M$^\dagger$instead of G$^2$M, further justification is required to substantiate the selection of the $\tau_q$ in Eq.3 (assuming it is the equation between lines 215 and 216) as the optimal threshold.

Essentially, the $\dagger$ version and the regular version differ only by the selection rule of $\tau_q$. The regular G$^2$M utilizes Eq. (3), and the G$^2$M$^\dagger$ considers the equation between lines 215 and 216. And according to [1], the two equations differ only when $\min_{t > 0} \lbrace t : \frac{1 + \vert \lbrace j : w_j \le -t \rbrace \vert}{\max(1, \vert \lbrace j : w_j \ge t \rbrace \vert)} \le q \rbrace > \min_{t > 0} \lbrace t : \sum_{j \in [p]} \mathbb{1} \lbrace j : w_j \ge t \rbrace < \frac{1}{q} \rbrace$, where $\mathbb{1}$ denotes the indicator function.

This means when $q$ is small or the fraction of non-nulls is small, then the first term usually finds a larger $\tau_q$, in favor of the FDR control. In this case, the power is pretty low, as when $\tau_q$ is large, there are only limited non-nulls selected. Instead of choosing $\tau_q$ with the first term (e.g., Eq. 3), the second term chooses a smaller $\tau_q$, which results in higher power. We will elaborate on the above in the revised manuscript near lines 215 and 216, where the new selection rule is introduced. In the revised manuscript, we will add an equation number to $\tau_q$ between lines 215 and 216.

QA3: It would be desirable to have at least one more real case study in section 3.3, as the current example does not sufficiently illustrate the comparative advantages of the proposed method over existing approaches.

We have prepared another real-world case study that focuses on a breast cancer dataset [2] that has 569 samples. It contains the binary label (i.e., cancer v.s. normal) as the dependent variable and 30 features as the independent variable. We have detailed the description of the data in the revised manuscript in Section 3.3.

Following a similar setting for the real case study, we investigated the literature and identified the features that have been proven to have high relevance to the cancer tissue. Using this as the ground truth, we report the “referenced/identified” features in the following table. In total, there are 22 relevant features.

It is clear that G$^2$M maintains the lowest false discovery rate, as there are only 3 missing. Meanwhile, it obtains 18 correctly identified features that are referenced, suggesting its good performance in real-world data.

The experimental setup is identical to the existing real case study; the only change is the dataset. And such information, including the results, is reflected in the revised manuscript. The selected features for each of the models are also explicitly displayed in the Appendix in the revised manuscript. We will add the results in the revised manuscript.

QA4: In section 3.1, the selection of $\rho$ needs some further justification.

We follow [3] to consider the $\rho=0.6$ setting. However, we also want to provide additional context on how the change of $\rho$ affects the feature selection performance. Hence, the $\rho=0.7$ setting is considered. We will add the above explanation in line 198 about this choice.

QA5: The simulation using the Copula model on page 6 needs additional clarification, given that the reader may not possess prior knowledge of this methodology.

We will include a new section in the Appendix in the revised manuscript to include the description of the Copulas. Essentially, Copulas are statistical tools designed to model and simulate complex dependencies among random variables, independently from the shapes of their marginal distributions. Unlike traditional multivariate models that assume linear or Gaussian relationships, copulas allow us to construct datasets where variables exhibit non-linear or asymmetric dependencies, better reflecting patterns seen in real-world data. In our study, we use two widely-studied copula families: the Clayton copula, which models strong lower-tail dependence (i.e., variables tend to move together when their values are low), and the Joe copula, which captures strong upper-tail dependence (i.e., variables tend to move together when their values are high). By specifying a copula parameter of 2, we control the strength of these dependencies in a consistent way across scenarios. For each simulated dataset, the individual variable distributions (marginals) are chosen to be either uniform (via identity transformation) or exponential (rate=1), allowing us to assess the robustness of our methods under different data distributions. This setup, implemented using the PyCop library, enables a comprehensive evaluation of how our proposed methods perform under diverse and realistic correlation structures.

QA6: Typos

Thanks for pointing out the typos. We will correct them in the revised manuscript.

References

[1] Ren, Z., & Barber, R. F. (2024). Derandomised knockoffs: leveraging e-values for false discovery rate control. Journal of the Royal Statistical Society Series B: Statistical Methodology, 86(1), 122-154.

[2] https://archive.ics.uci.edu/dataset/17/breast+cancer+wisconsin+diagnostic

[3] Sudarshan, M., Tansey, W., & Ranganath, R. (2020). Deep direct likelihood knockoffs. Advances in neural information processing systems, 33, 5036-5046.

QC5: The paper references paper 'Shen et al [32]' quite a lot, regarding problem set up, method, experiment setting, and result discussion. It would be good to add more discussion and clarifications on what's the novelty/originality of this paper compared with existing 'Shen et al [32]'

The FDR controlled feature selection problem is a widely studied problem in statistics and machine learning and its mathematical formulation is stated in Shen et al. and earlier works [1-8]. To address the problem, Shen et al. proposed a deep model-X knockoff approach. In this paper, we didn't build upon the method proposed in Shen et al. Instead, we extend the Gaussian mirror [8] and data splitting [2] methods, which do not need the generation of knockoff copies. The experiment setting in Shen et al. [32] gives a comprehensive coverage of existing experiment settings in FDR controlled feature selections [1-8]. Therefore, we follow the same experiment setting for a fair comparison with existing state-of-the-arts. We will clarify the difference and the connection to 'Shen et al [32]’ and other model-X-knockoff-based methods in the revised manuscript.

QC6: some discussion on the method complexity would be good and on how it compares with other methods in terms of computational complexity.

The complexity of the Gaussian Mirror method is $\mathcal{O}(np^3+p^4)$ (for $n>p$). Essentially it runs $p$ ordinary least square fit (OLS), each of which has $\mathcal{O}(np^2+p^3)$ complexity ($n>p$). The computational complexity of Algorithm 2 is $\mathcal{O}(np^3+p^4 + pki)$ (for $n>p$), where the $\mathcal{O}(pki)$ part is introduced by the k-means (following Lloyd's algorithm). $k$ is the number of clusters and $i$ stands for the number of iterations till convergence. In comparison, the data splitting runs two OLS, resulting in $\mathcal{O}(np^2+p^3)$ complexity. The knockoff framework, according to [10], requires at least $\mathcal{O}(p^{3.5})$ to solve for the knockoff variable in a semi-definite programming setting. Later on, it requires one OLS for $2p$ dimension, resulting in $\mathcal{O}(p^{3.5} + 4np^2+8p^3)$ complexity. Other deep-learning-based knockoff variables require additional deep learning model fitting to obtain the knockoff statistics, making the complexity analysis hard given the choice of optimizer and the model architecture, hence ignored. According to [11], the computational complexity for CRT, assuming $p=n$, is $\mathcal{O}(p^3 \log^2 p)$. We did not find any rigorous complexity analysis about HRT and dCRT, however, based on the results reported in Table 2 of the paper [12], we believe the complexity is between model-X knockoff and CRT. We will include this part in a separate section in the Appendix.

Thank you for the valuable feedback. Thank you for recognizing the strengths of our work, including a well-motivated problem with a clear rationale, a comprehensive discussion of related work, detailed theoretical proofs to support our method, and extensive experiments on both synthetic and real-world datasets against state-of-the-art approaches.

QC1: Novelty

We want to point out that we not only identified the limitation of the existing assumption for both Gaussian mirror and data splitting, but also derived the new test statistics given this relaxation. The benchmarking results demonstrated significant improvements, suggesting this relaxation on unit variance assumption is a strong cause of underperformance. Additionally, we also identified in the proof of Theorem 2.6 (i.e., section C.5 in the Appendix) that in order to boost power, one does not need to focus on jointly maximizing value of the non-null statistics.Instead, it is sufficient to consider treating each j-th test statistics separately. This observation is never explicitly pointed out by other papers, and is able to provide insights for future work that focuses on improving the power of feature selection algorithms with mirror test statistics. Overall, we believe all the above are valuable contributions of the paper. We will clarify and emphasize our contribution in the revised manuscript.

QC2: The experiment part can be enhanced with some significance test on the result compared with other methods.

Thank you for the suggestion. We want to point out that the false discovery rate (FDR) and the power are the gold standard for evaluating the performance of a feature selection algorithm, as indicated in [1-8]. Here, the power is defined as:

FDR is defined as:

And $\hat{S}_1 = \lbrace j: w_j \ge \tau_q\rbrace$.

QC3: It would be good to discuss more on how this normal distribution can be lifted/resolved.

Because we assume $\beta$’s are the coefficients for the linear model, so the fitting coefficients follow normal distributions [9]. Relaxing the normality assumption on the fitting coefficient essentially requires the relaxation on the linear model assumption. Instead of considering the relaxation on the normality assumption as the coefficients $\beta$’s (in the linear regression model) are simply parameters to quantify the relation between the input $X$ and the output $y$.

Our setting can be extended to a generalized linear model setting following the insights from [1]. However, the extension is non-trivial given that the change of the underlying model requires redefining the test statistics for G$^2$M. We will address this limitation and point out explicitly this direction as future work in section 4 of the revised manuscript.

QC4 : It would be good to add rationale on why k-means is selected in the estimation algorithm.

In Algorithm 2, the quantity of $\delta_j$, which is assumed to be given in Algorithm 1, needs to be estimated. In this practical setting where $\delta_j$ is not given, we treat the estimation of $\delta_j$ as an unsupervised learning problem given that we do not know the true values of $\delta_j$. And we choose the k-means given it has been widely considered as an efficient unsupervised algorithm. Specifically, we find the $k$ cluster centers from $\lbrace v_j\rbrace_{j=1}^{p}$. And the $k$ is chosen via the silhouette score (mentioned in line 164 in the manuscript). Then, we can give an estimate of $\delta_j$ by computing the deviation between $v_j$ and its cluster center as stated between line 10-12 in Algorithm 2. We will add the rationale in line 160 to link the estimation problem with k-means clustering.