Binh T Nguyen

Bertrand Thirion
email: bertrand.thirion@inria.fr

Sylvain Arlot
email: sylvain.arlot@universite-paris-saclay.fr

A Conditional Randomization Test for Sparse Logistic Regression in High-Dimension

Keywords:

Identifying the relevant variables for a classification model with correct confidence levels is a central but difficult task in high-dimension. Despite the core role of sparse logistic regression in statistics and machine learning, it still lacks a good solution for accurate inference in the regime where the number of features p is as large as or larger than the number of samples n. Here we tackle this problem by improving the Conditional Randomization Test (CRT). The original CRT algorithm shows promise as a way to output p-values while making few assumptions on the distribution of the test statistics. As it comes with a prohibitive computational cost even in mildly high-dimensional problems, faster solutions based on distillation have been proposed. Yet, they rely on unrealistic hypotheses and result in low-power solutions. To improve this, we propose CRT-logit, an algorithm that combines a variable-distillation step and a decorrelation step that takes into account the geometry of ℓ 1 -penalized logistic regression problem. We provide a theoretical analysis of this procedure, and demonstrate its effectiveness on simulations, along with experiments on large-scale brain-imaging and genomics datasets.

, with p as large as hundreds of thousands -compressible to thousands-but n stays at most few thousand. The family of methods

Introduction

Logistic regression is one of the most popular tools in modern applications of statistics and machine learning, partly due to its relative algorithmic simplicity. The method belongs to the class of generalized linear models that handle discrete outcomes, i.e. classification problems. Here, we focus on the binary classification problem, where one observation of the responses y ∈ {0, 1} and the data vectors x ∈ R p follows the relationship:

P(y = 1 | x) = g(x T β 0 ) = 1 1 + exp(-x T β 0 ) , ( 1
)
where g(x) = 1/(1 + exp(-x)) is the sigmoid function, and β 0 the vector of true regression coefficients. In the classical setting, in which the number of samples n is greater than the number of features p, an estimate β of the true signals β 0 can be obtained using maximum likelihood estimation (MLE). The asymptotic behaviour and derivation of the test statistic, confidence intervals and p-values of the MLE have been well studied, e.g. in [START_REF] Cox | Theoretical statistics[END_REF]. The availability of p-values for the test statistics makes it possible to rely on multiple hypothesis testing, where one wants to test which variables have a non-zero effect on the outcome, conditionally to the remaining variables, i.e.

(null) H j 0 : x j ⊥ y | x -j vs. (alternative) H j α :

x j ̸ ⊥ y | x -j ,
for each feature j ∈ [p]

def.

= {1, . . . , p} and x -j def.

= {x 1 , x 2 , . . . , x j-1 , x j+1 , . . . , x p }. Equivalently, under the setting in Eq. (1), we have:

(null) H j 0 : β 0 j = 0 vs. (alternative) H j α : β 0 j ̸ = 0. Unfortunately, this line of analysis cannot be applied to the high-dimensional regime, where p is larger than n, as argued in [START_REF] Sur | A modern maximum-likelihood theory for highdimensional logistic regression[END_REF]; [START_REF] Yadlowsky | SLOE: A Faster Method for Statistical Inference in High-Dimensional Logistic Regression[END_REF]; [START_REF] Zhao | The asymptotic distribution of the MLE in high-dimensional logistic models: Arbitrary covariance[END_REF]. These works show that in the regime where lim n,p→∞ n/p = κ, the MLE estimator exists only when κ > 2. However, we note that this type of analysis is done without the addition of ℓ 1 -regularization to the likelihood function, i.e. without using a penalized estimator to enforce sparsity. we consider is the Conditional Randomization Test (CRT) [START_REF] Candès | Panning for gold: 'model-x' knockoffs for high dimensional controlled variable selection[END_REF]. CRT relies on generating multiple noisy copies of original variables to output empirical pvalues in high-dimensional inference problems. However, prohibitive computational cost makes CRT impractical, as discussed at length in [START_REF] Candès | Panning for gold: 'model-x' knockoffs for high dimensional controlled variable selection[END_REF]; [START_REF] Tansey | The Holdout Randomization Test for Feature Selection in Black Box Models[END_REF]; [START_REF] Berrett | The conditional permutation test for independence while controlling for confounders[END_REF]; [START_REF] Liu | Fast and powerful conditional randomization testing via distillation[END_REF]. There have been several lines of research attempting to fix this problem, most notably the distilled Conditional Randomization Test (dCRT) [START_REF] Liu | Fast and powerful conditional randomization testing via distillation[END_REF]. This work introduced a distillation step as a replacement for the randomized sampling step to compute the importance statistics (see Section 2 for more details). It provides a way to output p-values for multiple types of regression and classification problems, assuming convergence to Gaussian of the test statistic in largesample regime. Yet, as shown in the left panel of Figure 1, the originally proposed dCRT test-statistic for logistic regression does not behave as well as intended. In particular, its null distribution deviates markedly from standard normal in high-dimension whenever n/p ≤ 1.

Contribution

We propose a correction for the dCRT, inspired by the decorrelation method presented in [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF]. The decorrelation step makes the nulldistribution of the test statistics much closer to standard normal, as shown on the right panel of Figure 1, and thus increases the statistical power of the method. We provide asymptotic analysis of this method, which shows that CRT-logit produces standard normal test-statistic in large-sample regime. In addition, we validate the high performance of CRT-logit on large-scale brain-imaging and genetics datasets, thus showing its usefulness in practical applications.

Related works

The closest cousin of the Conditional Randomization Test is Knockoff Filter [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF][START_REF] Candès | Panning for gold: 'model-x' knockoffs for high dimensional controlled variable selection[END_REF], a recent breakthrough in the False Discovery Rate (FDR) control literature. It relies on the creation of additional noisy features, called knockoffs, to calculate variable-importance statistics. Another extension of vanilla CRT is the Holdout Randomization Test (HRT) [START_REF] Tansey | The Holdout Randomization Test for Feature Selection in Black Box Models[END_REF]. While still requiring multiple sampling of noisy variables, HRT solves the computational issue of original CRT by doing heavy model fitting only once on one part of the dataset, and test statistics calculation on the other part, without refitting the model. However, this method relies on sample-splitting, hence inherently suffers from a loss of statistical power. A parallel line of work has introduced the Conditional Permutation Test (CPT) [START_REF] Berrett | The conditional permutation test for independence while controlling for confounders[END_REF], a non-parametric alternative to CRT that relies on a random shuffling mechanism applied to original variables, instead of multiple sampling of new variables. This potentially makes CPT more robust to model mis-specification. [START_REF] Yadlowsky | SLOE: A Faster Method for Statistical Inference in High-Dimensional Logistic Regression[END_REF] recently proposed a method called SLOE, which adapts the analysis of [START_REF] Zhao | The asymptotic distribution of the MLE in high-dimensional logistic models: Arbitrary covariance[END_REF], but in the regime different from what we are considering, where lim n,p→∞ n/p → κ ∈ (1, 2), and more importantly without sparsity-inducing penalty. On a separate note, we notice the similarity of dCRT [START_REF] Liu | Fast and powerful conditional randomization testing via distillation[END_REF] with debiased Lasso [START_REF] Javanmard | Confidence intervals and hypothesis testing for high-dimensional regression[END_REF][START_REF] Van De Geer | On asymptotically optimal confidence regions and tests for high-dimensional models[END_REF][START_REF] Zhang | Confidence intervals for low dimensional parameters in high dimensional linear models[END_REF]. This line of work proposed a debiasing formula for the estimator, which makes the asymptotic distribution of ( βLASSOβ 0 ) standard normal, so that one can compute the test statistic and p-value associated with each variable.

Background

Notation We denote matrices, vectors, scalars and sets by bold uppercase, bold lowercase, script lowercase , and calligraphic letters, respectively, e.g. X, x, x, and X . The i-th row of a matrix X will be denoted X i, * , the j-th column X * ,j and the (i, j)-th element X i,j . For any natural number p, we denote the set [p] def.

= {1, . . . , p}. For each x ∈ R p and j ∈ [p], we denote x -j def.

= {x 1 , x 2 , . . . , x j-1 , x j+1 , . . . , x p } a p-1 dimensional observation after removing the j-th variable. Correspondingly, X -j is the data matrix X ∈ R n×p with column X * ,j removed. The cumulative distribution function (CDF) of the standard Gaussian distribution will be denoted Φ(•). The indicator function of a random event A will be denoted 1 A . For any two positive sequences x n and y n , we write x n ≍ y n if cy n ≤ x ≤ Cy n for all n, for some positive constants c and C. For a vector x, ∥x∥ p denotes its ℓ p norm. For a function f : R p → R, ∇ j f denotes its gradient w.r.t. the j-th variable, for j ∈ [p].

Problem setting

We consider exclusively binary classification, where the response vector is denoted y ∈ {0, 1} n and the data matrix X ∈ R n×p consists of n p-dimensional samples. Throughout the paper, we assume the data {X i, * } n i=1 are i.i.d. and follow a distribution with zero mean and population covariance matrix Σ. Moreover, we assume that X i, * and y i follow the logistic relationship in Eq. ( 1). We denote the support set S def.

= {j ∈ [p] : β 0 j ̸ = 0} and assume that it is sparse, i.e. card(S) = s * ≪ p, where card denotes the cardinality of a set. Furthermore, Ŝ def. = {j ∈ [p] : βj ̸ = 0} indicates an estimation of S, where βj is an estimate of the true signal β 0 j . We try to obtain it through a ℓ 1 -penalized logistic estimator:

βλ = argmin β∈R p ℓ(β) + λ∥β∥ 1 , with ℓ(β) = - 1 n n i=1 {y i (X i, * β) -log [1 + exp(X i, * β)]} .
(2)

We denote I def.

= E β 0 [∇ 2 ℓ(β 0 )] the Fisher information matrix, and I j|-j the partial Fisher information, defined by

I j|-j def. = E[∇ 2 jj ℓ(β 0 )-[∇ 2 j,-j ℓ(β 0 )] ⊤ [∇ 2 -j,-j ℓ(β 0 )] -1 ∇ 2 -j,j ℓ(β 0 )] = I jj -I j,-j I -1
-j,-j I -j,j , where I j,-j is the row-vector made with the jth-row and the columns corresponding to β -j , I -j,-j the sub-matrix of I made with the rows and columns corresponding to β -j . This quantity, defined following Cox and Hinkley (1979, pp. 323), plays an important role in our proposed method, detailed in Section 3.

Statistical control with False Discovery Rate

To quantify statistical errors, we consider the False Discovery Rate, introduced in [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF]. Given an estimate of the support Ŝ, the false discovery proportion (FDP) is the ratio of the number of selected features that do not belong to the true support S, divided by the total number of selected features. The False Discovery Rate is the expectation of the FDP:

FDP( Ŝ) = card({j : j ∈ Ŝ, j / ∈ S}) card( Ŝ) ∨ 1 and FDR( Ŝ) = E[FDP( Ŝ)].

Conditional Randomization Test (CRT) and Distillation CRT (dCRT)

The concept of Conditional Randomization Test was originally proposed in the model-X knockoff paper [START_REF] Candès | Panning for gold: 'model-x' knockoffs for high dimensional controlled variable selection[END_REF] as a way to output valid empirical p-values using knockoff variables. The principle of the knockoff filter is first to sample noisy copies X * ,j of variable X * ,j , given a known sampling mechanism P j |-j . One advantage of the knockoff filter is that no specific assumption is placed on the distribution of the inferred test statistic. However, this means that, in general, there is no mechanism to derive p-values from the knockoff statistic. This motivates the introduction of CRT, which requires running high-dimensional inference for each variable j B times. However, the computation cost of CRT is prohibitive when p grows large: assuming that we use the Lasso program with coordinate descent to compute T CRT j , its runtime would be O(Bp 4 ) (Hastie et al., 2009, pp. 93). Moreover, CRT requires decently large B to make the empirical distribution of the p-values smooth enough. Reducing the computational cost of CRT is the main motivation of several works [START_REF] Berrett | The conditional permutation test for independence while controlling for confounders[END_REF][START_REF] Liu | Fast and powerful conditional randomization testing via distillation[END_REF][START_REF] Tansey | The Holdout Randomization Test for Feature Selection in Black Box Models[END_REF]. One of them is the introduction of distillation-CRT (dCRT) by [START_REF] Liu | Fast and powerful conditional randomization testing via distillation[END_REF]. The main appeal of this method is that it can output p-values analytically, therefore bypassing the multiple knockoffs sampling steps used to infer on each variable, and leads to a reasonable reduction of the computation cost.

Distillation operation

The key addition of dCRT is the distillation operation: for each variable j, we want to distill all the conditional information of the remaining variables X -j to x j and to y via least-squares minimization with ℓ 1 -regularization to enforce sparsity. For each variable j, we first solve the lasso problem by regressing X * ,j on X -j , βd X * ,j = argmin

β∈R p-1 1 2 ∥X * ,j -X -j β∥ 2 2 + λ dx ∥β∥ 1 . ( 3
)
For distillation of variable j and the binary response y with logistic relationship, [START_REF] Liu | Fast and powerful conditional randomization testing via distillation[END_REF] briefly suggested to solve a penalized estimation problem, similar to Eq. ( 2):

βdy,j = argmin

β∈R p-1 - 1 n n i=1 y i (X ⊤ i,-j β) -log 1 + exp(X T i,-j β) + λ∥β∥ 1 . (4)
Then, a test statistic is calculated for each j = 1, . . . , p:

T j = √ n ⟨y -X -j βdy,j , X * ,j -X -j βd X * ,j ⟩ ∥y -X -j βdy,j ∥ 2 ∥x j -X -j βd X * ,j ∥ 2 , ( 5
)
which, under the null hypothesis, and more importantly, assuming linear relationship between X i, * and y, follows standard normal distribution asymptotically, conditional to y and X -j . It then follows that we can output a p-value for each variable j by taking

p j = 2 [1 -Φ (T j )].
However, the formulation of test statistics in Eq. ( 5) is not truly satisfactory in the setting of sparse logistic regression. More specifically, both the calculation of regression residuals y -X -j βdy,j and test statistics T j do not take into account the non-linear relationship between X and the binary response y. The first row of Figure 6 plots the qq-plot of the test statistics T j for logistic regression, which shows that even in the classical regime where n > p, its distribution is far from standard normal.

Decorrelating Test-Statistics for High-Dimensional Logistic Regression

As we have elaborated, the formulation of dCRT is not well-suited for problems other than penalized least-squares regression. We therefore propose an adaptation of dCRT in the case of logistic regression, inspired by the classical work of [START_REF] Cox | Theoretical statistics[END_REF] and by [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF]. First, note that when testing H j 0 : β 0 j = 0 under the case where n > p, we have the classical Rao's test statistic, defined by

T Rao j = √ n Î-1/2 j|-j ∇ j ℓ(0, β-j ) , ( 6
)
where

∇ j ℓ(0, β-j ) def. = ∇ βj ℓ(β j , β-j ) βj =0
is the Fisher score. Here β-j def.

= argmin

β-j ∈R p-1 ℓ(β j , β -j )
is the constrained maximum-likelihood estimator of β -j with fixed β j , and Îj|-j is a consistent estimator of the partial Fisher information I j|-j . The appearance of the term Î-1/2 j|-j is due to the fact that under the null hypothesis H j 0 , we have, by Cox and Hinkley (1979, Chapter 9), and by [START_REF] Rao | Large sample tests of statistical hypotheses concerning several parameters with applications to problems of estimation[END_REF],

√ n∇ j ℓ(0, β-j ) (d) ----→ n→∞ N (0, I j|-j ) ,
which makes the asymptotic distribution of T Rao j standard normal. However, in the highdimension case, where n < p, we do not reach this convergence in distribution. To see this, consider the Taylor expansion of the Fisher score of variable j around any given estimator β -j of the true β 0 -j :

∇ j ℓ(0, β -j ) = ∇ j ℓ(0, β 0 -j ) + ∇ 2 j,-j ℓ(0, β 0 -j )( β -j -β 0 -j ) + O ( β -j -β 0 -j ) 2 (7)
On the right-hand side, the first term converges weakly to a normal distribution due to the Central Limit Theorem, the remainder term becomes negligible using ℓ 1 regularization to induce sparsity, but the second term does not, due to estimation bias and sparsity effect of β -j [START_REF] Fu | Asymptotics for lasso-type estimators[END_REF].

Adapting distillation operation for sparse logistic regression Fortunately, Eq. ( 7) suggests that for each variable j, we can debias the Fisher score by correcting the impact of other terms. In particular, for each variable j, we replace the Fisher score by

∇ j ℓ(0, β -j ) -I j,-j I -1 -j,-j ∇ -j ℓ(0, β -j ) . ( 8
)
The inversion of the large matrix I -j,-j ∈ R (p-1)×(p-1) is computationally prohibitive, but we can estimate the term I j,-j I -1 -j,-j straightforwardly by solving ŵj = argmin

w∈R p-1 1 2n n i=1 ∇ 2 j,-j ℓ i ( β) -w T ∇ 2 -j,-j ℓ i ( β) 2 + λ∥w∥ 1 , (9)
for each variable j, where β is given with Eq. ( 2). Moreover, since we have the closedform of the derivatives of the logistic loss ℓ( β), a simple derivation from Eq. ( 9) suggests the following x j -distillation operation, adapted for logistic regression:

βd X * ,j = argmin β∈R p-1 1 n n i=1 g ′′ (X i, * β)(X i,j -X i,-j β) 2 + λ dx ∥β∥ 1 , ( 10
)
where the extra term (second-order derivative of the sigmoid function) g ′′ (X i, * β) = exp (Xi, * β)

[1+exp (Xi, * β)] 2 appears from differentiating twice the loss function ℓ( β), and β = βλ is defined in Eq. ( 2). On the other hand, we can obtain βdy,j j from β by simply omitting the j-th coefficient from it, i.e. βdy,j def.

= ( β1 , β2 , . . . , βj-1 , βj+1 , . . . , βp ) .

Finally, the equation for decorrelated test score, adapted from both Eq. ( 5) and ( 6), reads

T decorr j = - 1 √ n Î-1/2 j|-j n i=1 y i -g(X i,-j βdy,j ) X i,j -X i,-j βd X * ,j , ( 11
)
where the formula for empirical partial Fisher information is Îj|-j

= n -1 n i=1 g ′′ (X i, * β)(X i,j - X i,-j
βd X * ,j ) X i,j . A summary of the full procedure, which we call CRT-logit, can be found in Algorithm 1. Notice that the runtime of CRT-logit is the same as dCRT, which means in general slower than KO and HRT. To speedup inference time, we introduce a variable-screening step that eliminates potentially unimportant variables before distillation, similar to dCRT. We provide empirical benchmark of runtime of each method in Section 4.5.

Setting ℓ 1 -regularization parameter λ and λ dx In general, we advise to use crossvalidation for obtaining βλ in Eq. ( 2) and for X * ,j -distillation operator, as defined by Eq. ( 10). This is inline with the theoretical argument for dCRT [START_REF] Liu | Fast and powerful conditional randomization testing via distillation[END_REF], Lemma 1 and Theorem 3). However, we also observe empirically that choosing the ℓ 1regularization parameters of the distillation step can strongly affect how variables are selected by CRT-logit. We provide more details Appendix C, and leave further theoretical investigations of this phenomenon for future work.

Algorithm 1: CRT-logit

1 INPUT design matrix X ∈ R n×p , reponses y ∈ R n 2 OUTPUT vector of p-values {pj} p j=1 ; 3 βλ ← solve_sparse_logistic_cv(X, y) // Using Eq. (2) 4 ŜSCREENING ← {j : j ∈ [p], βj ̸ = 0} 5 for j ∈ ŜSCREENING do 6 βd X * ,j ← solve_scaled_lasso_cv(X * ,j, X * ,-j) //

Asymptotic analysis of the Decorrelated Test Statistic

We now provide theoretical analysis of CRT-logit in large-sample regime. All the proofs can be found in Appendix A. Without writing it explicitly, in our analysis, we consider p = p(n), and the following assumption. Assumption 3.1 (Regularity conditions). Assume that (A1) λ min (I) ≥ κ 2 for some constant κ > 0. (A2) Sparsity of β 0 and w 0,j , with w 0,j the ground truth weights for the distillation of x j in Eq. (10): |S| = s * and ∥w 0,j ∥ 0 = s ′ with s * = o n 1/2 / log(p) and

s ′ = o n 1/2 / log(p) . (A3) For all i ∈ [n], X i, * and (-y i + g ′ (X i, * β)) are sub-exponential random variables,
and |X i,-j w 0,j | ≤ K almost surely, for some constant K.

We then have the following result, that states that the asymptotic distribution of the decorrelated test scores is standard normal.

Theorem 3.1. Let j ∈ [p], and let T decorr j be defined as in Eq. (11), with λ ≍ λ dx ≍ n -1 log(p). Then, if Assumption 3.1 holds true, under the null hypothesis H j 0 :

β 0 j = 0, we have ∀t ∈ R , lim n→∞ |P β 0 (T decorr j ≤ t) -Φ(t)| = 0 ,
where Φ(•) is the CDF of the standard Gaussian distribution. Moreover, for each j ∈ [p],

if we define p j def.

= 2 [1 -Φ (T j )] , i.e. p j is the output of Algorithm 1, then, under the null hypothesis H j 0 : β 0 j = 0, we have

lim sup n→∞ P β 0 ( p j ≤ t) ≤ t for all t ∈ [0, 1] ,
that is, the p-values output by Algorithm 1 are valid asymptotically.

FDR control with CRT-logit

We now state the second main result, which establishes that the FDR of the test is controlled when using Benjamini-Yekutieli procedure [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] with the p-values output from Algorithm 1, assuming that the number of tests p is fixed.

Theorem 3.2. Under Assumptions 3.1 and logistic model defined in Eq (1), with λ ≍ λ dx ≍ n -1 log(p), assume n -1/2 (s ′ ∨ s * ) log(p) = o(1), and assume the number of tests p is fixed. Let α ∈ (0, 1) and S BY-CRT be given by applying following the Benjamini-Yekutieli FDR-controlling procedure to the CRT-logit p-values { p j } j∈ [p] , output from Algo.1. Then, we have

lim sup n→∞ E card( S BY-CRT ∩ S c ) card( S BY-CRT ) ∨ 1 ≤ α .
Remark 3.1. Assumption 3.1 is also assumed in [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF]; van [START_REF] Van De Geer | On asymptotically optimal confidence regions and tests for high-dimensional models[END_REF], which also provide a detailed discussion of this regularity assumption in generalized linear models. This assumption, in turn, is built on the regularity assumption in the classic work of Cox and Hinkley (1979, Chapter 9) to establish asymptotic normality of Rao's test statistic. Theorem 3.1 is an adaptation of Ning and Liu (2017, Theorem 3.1), specialized for the case of sparse logistic regression and the p-values output from CRT-logit.

Empirical Results

We provide benchmarks of the proposed CRT-logit algorithm along with most other methods mentioned in the introduction, in particular: model-X Knockoff (KO) [START_REF] Candès | Panning for gold: 'model-x' knockoffs for high dimensional controlled variable selection[END_REF], Debiased Lasso (dLasso) [START_REF] Zhang | Confidence intervals for low dimensional parameters in high dimensional linear models[END_REF][START_REF] Javanmard | Confidence intervals and hypothesis testing for high-dimensional regression[END_REF], original CRT with 1000 samplings [START_REF] Candès | Panning for gold: 'model-x' knockoffs for high dimensional controlled variable selection[END_REF], Holdout Randomization Test with 5000 samplings [START_REF] Tansey | The Holdout Randomization Test for Feature Selection in Black Box Models[END_REF], and Lasso-Distillation CRT (dCRT) [START_REF] Liu | Fast and powerful conditional randomization testing via distillation[END_REF]. We did not include SLOE [START_REF] Yadlowsky | SLOE: A Faster Method for Statistical Inference in High-Dimensional Logistic Regression[END_REF] and CPT [START_REF] Berrett | The conditional permutation test for independence while controlling for confounders[END_REF], as the provided open-source implementation are particularly unstable and do not fit in the sparse-regression setting (for SLOE), or implementation are not available (for CPT).

Remark 4.1. As a slight caveat, in the simulated and semi-realistic experiment sections (Sections 4.1,4.2 and 4.3), we introduce an additional noise term to the logistic relationship of Eq. (1):

P(y i = 1 | x i ) = g(x T i β 0 + σξ i ) , ( 12
)
where ξ i ∼ N (0, 1) is a Gaussian noise and σ > 0 the noise magnitude. The formula in Eq. ( 12) has been used in previous works, e.g. [START_REF] Bzdok | Semi-Supervised Factored Logistic Regression for High-Dimensional Neuroimaging Data[END_REF]. There is a clear justification to this: in most of the applications we consider, data are collected with measurement errors. In the case of brain-imaging, for example, recording the brain signal of the human subjects by scanners often includes noise caused either from the machine, or from the movement of the subjects, as elaborated by [START_REF] Lindquist | The Statistical Analysis of fMRI Data[END_REF]. Moreover, in general, this setting corresponds to a model mis-specification, which the CRT-logit is robust to under Assumption 3.1, following the arguments as in Ning and Liu (2017, Section 5).

Remark 4.2. We use Benjamini-Hochberg step-up procedure [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF] to control FDR with the p-values in all the empirical experiments in Section 4.2 and App. 4.4, as we observe empirically the FDR is usually controlled with this procedure, without compromising power with the conservative BY bound.

Effectiveness of the decorrelation step

To show how decorrelating the test statistics helps, we set up a simulation with matrix X of p = 400 features and vary the number of samples n ∈ {200, 400, 800, 4000}. The binary response vector y is created following Eq. ( 12), and the design matrix X is sampled from a multivariate normal distribution with zero mean, while the covariance matrix Σ ∈ R p×p is a symmetric Toeplitz matrix, where the parameter ρ ∈ (0, 1) controls the correlation between neighboring features: correlation decreases quickly when the distance between feature indices increases. The true signal β 0 is picked with a sparsity parameter κ = s * /p that controls the proportion of non-zero elements with magnitude 2.0, i.e. β j = 2.0 for all j ∈ S. For the specific purpose of this experiment, non-zero indices of S are kept fixed. The noise ξ is i.i.d. normal N (0, Id n ) with magnitude σ = ∥Xβ 0 ∥ 2 /( √ n SNR), controlled by the SNR parameter. In short, the three main parameters controlling this simulation are correlation ρ, sparsity degree κ and signal-tonoise ratio SNR. We generate randomly 1000 datasets, and run dCRT and CRT-logit algorithm to obtain one sample of test statistics {T j } p j=1 and {T decorr j } p j=1 . Then, we pick 1000 samples of one null test statistic T j and T decorr j , defined in Eq. ( 5) and ( 11), respectively, and plot the qq-plot of their empirical quantile versus the standard normal quantile. From the results in Figure 6, we observe that the empirical null distribution of the test statistic is much closer to a standard normal when adding the decorrelation step.

In particular, when the sample size n increases to 400, the decorrelated test statistic has empirical quantiles almost inline with the theoretical quantiles of the standard normal distribution, while dCRT test score strays away from the 45-degree line. Again, we emphasize that the normality of T j is crucial for the p-values calculation. This outlines the importance of the decorrelating step on T j .

High-dimensional scenario with varying simulation parameters

To see how each algorithm performs under different settings, we follow the same simulation scenario as in Sec. 4.1, but vary each of the three simulation parameters, while keeping the others unchanged at default value of SNR = 2.0, ρ = 0.5, κ = 0.04. We target a control of FDR at level 0.1, using Benjamini-Hochberg procedure. Results in Figure 3 show that CRT-logit is the most powerful method while still controlling the FDR. Moreover, in the presence of higher correlations between nearby variables (ρ > 0.6), other methods suffer a drop in average power, but this is not as severe for CRT-logit. The original CRT, in general, is conservative. We believe that this is due to using only B = 500 samplings to generate empirical p-values for the two methods, due to prohibitive average runtime of the algorithm with larger B (which we provide in Section 4.5). For HRT, the conservativeness is expected, due to the usage of only half of the sample for test-statistics calculation -even though the number of samplings is bigger than original CRT (B = 5000). We note that, perhaps surprisingly, the debiased lasso (cdlasso) is the most conservative. It controls FDR well in all settings. This might be due to the fact that dlasso also relies on the choice of the ℓ 1 -regularization λ in the nodewise Lasso operation, similar to the X * ,j -distillation of dCRT, as noted in Section 1. What makes the difference is that instead of using cross-validation for setting λ for each variable j, a fixed value of λ = 10 -2 λ max is used in the implementation of dlasso. We strongly suspect this fixed value is not optimal, which makes the procedure powerless.

Application: large-scale analysis on brain-imaging dataset

Description The Human Connectome Project dataset (HCP) is a collection of brain imaging data on healthy young adult subjects with age ranging from 22 to 35. More specifically, the input X is a set of 2mm statistical maps of 400 subjects across 8 cognitive tasks. These are called z-maps, as the data follow a standard normal distribution under the null hypothesis. Each task in turn features 2 different contrasts, which effectively form binary responses y ∈ {0, 1} n . We propose to fit y through distributed brain signals and identify relevant brain locations. The setting is high-dimensional with n = 800 samples, corresponding to 400 subjects, while the total number of variables is p ≈ 200, 000 brain voxels. Following [START_REF] Nguyen | Ecko: Ensemble of clustered knockoffs for robust multivariate inference on fmri data[END_REF]; [START_REF] Chevalier | Decoding with confidence: Statistical control on decoder maps[END_REF], we use a hierarchical clustering scheme to group the variables into C = 1000 spatially connected clusters. We provide details of the pre-processing step in Appendix F.

Creating semi-realistic ground-truth and response labels Since there is no ground truth for this dataset, we create synthetic true signals by fitting the data X and response y with an ℓ 1 -penalized logistic classifier. In other words, the estimator βlogreg will serve as true regression coefficients for each task. Then, to avoid bias in simulating label ŷ, the z-maps matrix X of one task are used in conjunction with the discriminative pattern map βlogreg from the next task in the following order: relational, gambling, emotion, social. For instance, we use βlogreg of gambling with z-maps data matrix of relational, i.e. for all i = 1, . . . , n, given

x i,relational , ŷi ∼ Bern g(x ⊤ i,relational βlogreg gambling + σξ i ) , ( 13
)
where Bern(a) is a Bernoulli probability mass function that takes a value 1 with probability a, σ is a noise magnitude and ξ i is a standard normal noise. Finally, we apply all inference algorithms on the semi-synthetic data (X, ŷ), and we evaluate their performance using the ground-truth βlogreg . This simulation setting is similar to [START_REF] Chevalier | Decoding with confidence: Statistical control on decoder maps[END_REF], except that here we consider a classification and not a regression problem. It allows us to calculate the False Discovery Rate and average power with multiple runs of the inference procedure (across tasks).

Remark 4.3. The i.i.d. assumption is formally violated in this experiment, where for each subject we analyze two sample image that are not independent. Yet, this remains a short-range correlation structure, and is thus not a strong challenge to the i.i.d. assumption.

Results

Results in Figure 4 show that CRT-logit achieves a better recovery compared to KO or original CRT/dCRT/HRT, which results in higher statistical power. This gain comes with a good control of the FDR under desired level α = 0.1. On a related note, the only analysis where dCRT makes more discoveries than CRT-Logit is in emotion task, but at the cost of failing to control FDR at nominal level.

Application: genome-wide association study with Human Brain Cancer Dataset

Description The last in our benchmark is a Genome-wide Association Study (GWAS) on the The Cancer Genome Atlas (TCGA) dataset [START_REF] Weinstein | The cancer genome atlas pan-cancer analysis project[END_REF][START_REF] Vasaikar | LinkedOmics: analyzing multi-omics data within and across 32 cancer types[END_REF]. We choose to analyze the Glioma cohort, which consists of n = 1026 patients across a wide age range, diagnosed with this type of brain tumor, with a total of p = 24776 genes in the data matrix, recorded as copy number variations (CNVs) at the gene level in log ratio format. As with the brain-imaging inference in Section 4.3, we use clustering to reduce the dimension to C = 1000 clusters. For the response, a long-term survivor (LTS) is defined as a patient who survived more than five years after diagnosis and would be labeled y = 0, and any patient who died within five years would be a shortterm survivor (STS), labeled y = 1. The objective is to identify significant genes that contribute to classification of the LTS/STS status. Similar to the Human Connectome Project dataset, there is no real ground-truth for the TCGA Glioma. However, we have the list of mutations and the frequency of those detected in the diagnosed patients. We therefore select the 1000 most frequent gene mutations that appeared in this list, i.e. the ground truth list consists of 1000 genes (variables).

Result

The result from Table 1 shows that CRT-logit finds the largest number of genes. Moreover, most of selected genes in this table are detected in the list of mutated genes found on recorded patients. Some genes are detected by all the benchmarked methods, most prominently SPEN, which is found on over 10 % of patients in the cohort. [START_REF] Vasaikar | LinkedOmics: analyzing multi-omics data within and across 32 cancer types[END_REF]. Furthermore, this gene is known to be associated not only with brain cancer, but also with other types of cancer in The Human Protein Atlas project [START_REF] Légaré | The estrogen receptor cofactor spen functions as a tumor suppressor and candidate biomarker of drug responsiveness in hormone-dependent breast cancers[END_REF]. Note that, in the absence of a ground-truth, this does not guarantee all genes found are associated with glioma, but this experiment demonstrates the capability of CRT-logit in GWAS studies. Besides statistical performance, it is equally important to assess the computational cost of inference procedures. The average runtime in Table 2 from the two experiments shows that the original CRT is not suitable for large-scale inference: it is over 2000 times slower than the fastest method (Knockoff Filter), and over 150 times slower than dCRT/CRT-logit. The empirical runtime also confirms the effectiveness of the screening step before doing distillation/decorrelation of the test-statistics: the step makes CRTlogit and dCRT 20 times faster than without. On a related note, although in theory Debiased Lasso, dCRT and CRT-logit (both without screening) share the same runtime complexity, the latter two are slower due to the use of cross-validation to estimate the sparsity hyperparameter λ and λ dx (detailed in Section 3).

Methods

Average runtime of benchmarked methods

Discussion

We proposed an adaptation of the Conditional Randomization Test (CRT) for sparse logistic regression in the high-dimensional regime. A major improvement of CRT-logit, our proposed algorithm, compared to original CRT, comes from the decorrelation of test statistics to make their distribution closer to standard normal. Indeed, results from synthetic experiments in Figure 6 show that in high-dimension (when 0.5 ≤ n/p ≤ 1.0), the empirical null distribution of CRT-logit's test statistic T decorr is much more similar to a standard normal compared to the original CRT test statistic. Moreover, empirical benchmarks in Section 4 demonstrate that CRT-logit performs better than related statistical inference methods, such as the Debiased Lasso or Model-X Knockoffs. In particular, CRT-logit is the most powerful method in our synthetic experiment with high-dimensional datasets in Section 4.2, while still keeping FDR controlled under predefined level α = 0.1.

We note that there exists some limitations to CRT-logit. The computational cost of CRT-logit, while lower than vanilla CRT, is still larger than alternative methods such as Knockoff Filter and Holdout Randomization Test. Moreover, tuning the ℓ 1 -regularization λ dx parameter by cross-validation, as is often done, can further increase the computational cost of CRT-logit (and dCRT).

Despite this, our empirical benchmarks on both simulated and real data show real promises of CRT-logit. Henceforth, we believe CRT-logit is competitive for practical settings that involve structured data, such as brain-imaging and genomics applications.

where we use triangle inequality with the last step. By Taylor expansion, and from Lemma A.4, we have

A 1 ≤ n 1/2 |v * ⊤ ∇ 2 ℓ(β 0 )( βnull -β 0 )| ≤ n 1/2 ∥ βnull -β 0 ∥ 1 ∥v * ⊤ ∇ 2 ℓ(β 0 )∥ ∞ ≾ s * log(p) √ n
where the second inequality is by Holdër inequality, and the last inequality us due to Lemma A.1 and A.2. Similarly, we can bound A 2 , by using Lemma A.3 and Lemma A.4

A 2 ≤ n 1/2 |(v -v * ) ⊤ ∇ℓ(β 0 )| ≤ n 1/2 ∥v -v * ∥ 1 ∥∇ℓ(β 0 )∥ ∞ ≾ (s * ∨ s ′ ) log p n
This implies that,

n 1/2 |v ⊤ ∇( βnull ) -v * ⊤ ∇ℓ(β 0 )| ≾ n -1/2 (s * ∨ s ′ ) log(p) . ( 15
)
The remaining part of the proof is to bound Îj|-j -I j|-j , where, by definition

I j|-j = E g ′′ (X i, * β 0 ) X i,j -X i,-j w 0,j X i,j
Evaluating the difference between Îj|-j and I j|-j gives Îj|-j -

I j|-j = 1 n n i=1 g ′′ (X i, * β) X i,j -X i,-j βd X * ,j X i,j -E g ′′ (X i, * β 0 ) X i,j -X i,-j w 0,j X i,j = 1 n n i=1 g ′′ (X i, * β)X 2 i,j -E g ′′ (X i, * β 0 )X 2 i,j + 1 n n i=1 g ′′ (X i, * β)X i,-j βd X * ,j X i,j -E g ′′ (X i, * β 0 )X i,-j w 0,j X i,j ≤ 1 n n i=1 g ′′ (X i, * β)X 2 i,j -E g ′′ (X i, * β 0 )X 2 i,j C + 1 n n i=1 g ′′ (X i, * β)X i,-j ( βd X * ,j -w 0,j ) X i,j B1 + 1 n n i=1 [g ′′ (X i, * β 0 ) -g ′′ (X i, * β)]X i,-j w 0,j X i,j B2 + 1 n n i=1 g ′′ (X i, * β 0 )X i,-j w 0,j X i,j -E g ′′ (X i, * β 0 )X i,-j w 0,j X i,j B3 ,
where the last step follows triangle inequality.

We have, by Cauchy-Schwartz inequality, by Lemma A.3, and by the fact that g ′′ (x) ∈ (0, 1) for every x ∈ R; and X i,-j , X i,j is sub-exponential by Assumption 3.1:

B 1 ≤ 1 n n i=1 g ′′ (X ⊤ i, * β)(( βd X * ,j -w 0,j ) ⊤ X i,-j ) 2 1 n n i=1 g ′′ (X ⊤ i, * β)X 2 i,j ≾ (s * ∨ s ′ ) log(p) n .
Similarly, to bound B 2 , we have, again by Cauchy-Schwartz inequality,

B 2 ≤ 1 n n i=1 [g ′′ (X i, * β 0 ) -g ′′ (X i, * β)] 2 (X i,-j w 0,j X i,j ) 2 ≤ 1 n n i=1 [g ′′ (X i, * β 0 )X i, * ( β -β 0 )] 2 (X i,-j w 0,j X i,j ) 2 ,
where the second inequality comes from using the self-concordance property of the sigmoid function (discussed at length in [START_REF] Bach | Self-concordant analysis for logistic regression[END_REF] This implies

| Îj|-j -I j|-j | ≾ (s * ∨ s ′ ) log p n . ( 16
)
Putting Equation ( 15) and ( 16) together, we have, under null hypothesis,

T decorr j D -→ n 1/2 I -1/2 j|-j v * ⊤ ∇ℓ(β 0 ) def. = T * j ,
with convergence rate O(n -1/2 ). Finally, by noting that we can decompose ∇ℓ(β 0 ) = 1 n n i=1 ∇ℓ i (β 0 ), and each ℓ i (β 0 ) has bounded first, second, and third moment, a direct application of Berry-Esseen theorem give convergence in distribution of T * j to a standard normal law, with rate O(n -1/2 ).

We also arrive at the second conclusion of Theorem 3.1 by noting that it is a straightforward by-product of the result on normality of the distribution of decorrelated test score under null hypothesis, based on the formula for the p-values of CRT-logit algorithm.

Proof of Theorem 3.2. The proof of this theorem is a straightforward adaptation from [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]. For shorter notation, we denote S = i∈S c j j ′ =1 p k=1 p i,j ′ ,k for all j ∈ {1, . . . , p}, and remark that p i,j ′ ,k = 0 if j ′ > k, by definition of S BY-CRT . We then have

A = F (1) + p j=2 1 j F (j) -F (j -1) = p-1 j=1 1 j - 1 j + 1 F (j) + F (p) p .
This leads to

E | S ∩ S c | | S| ∨ 1 ≤ p-1 j=1 1 j - 1 j + 1 F (j) + F (p) p (19)
By the definition of p i,j,k in Eq. ( 18), we have (Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995)). Let α ∈ (0, 1) be the predefined FDR control level. Let p1 , . . . , pm be output p-values from inference algorithm, e.g. Algorithm 1. We reorder them ascendingly, denoted by p(1) ≤ We omit Holdout Randomization Test (HRT) as the method does not select any brain region. For dlasso, dCRT and CRT-logit, we plot the test-statistics; for KO the sign of selected coefficients, and for CRT the -log 10 of the empirical p-values.

F (j) =

Figure 1 :

 1 Figure 1: QQ-Plot for 1000 samples of test-statistic of a null index for logistic regression, with simulated data, n = 200, p = 400. Left: Statistics obtained from running Distilled-CRT, and Right: from our proposed algorithm. The empirical distribution of dCRT null-statistic strays far from theoretical distribution, which is standard normal, while empirical distribution of CRT-logit's null test score is much closer.

Figure 2 :

 2 Figure 2: QQ-Plot for one null CRT statistic for logistic regression, with varying number of samples and a fixed number of variables p = 400. The theoretical quantiles are obtained from a standard Gaussian distribution. The decorrelation step makes the empirical null distribution of the null statistics much closer to standard Gaussian. Parameters: SNR = 3.0, ρ = 0.4, sparsity = 0.06. Upper row: Distilled-CRT statistic defined by Eq. (5). Bottom row: CRT-logit, with decorreleated test score defined by Eq. (11) (ours).

Figure 3 :

 3 Figure 3: FDR/Average Power of 100 runs of simulations across varying parameters in high-dimensional settings. Default parameter: n = 400, p = 600, SNR = 2.0, ρ = 0.5, κ = 0.04. FDR is controlled at level α = 0.1. Methods: Debiased Lasso (dlasso), model-X Knockoff (KO-logit), original CRT (CRT), HRT (HRT), dCRT (dCRT), and our version of CRT (dark green line -CRT-logit).

Figure 4 :

 4 Figure 4: FDR/Average Power of 50 runs of semi-realistic experiments on four tasks of Human Connectome Project dataset. Parameters: n = 800 (taken from 400 subjects), SNR = 2.0. Methods (clustering versions): Debiased Lasso (cdlasso), model-X Knockoff (cKO-logit), original CRT (CRT), HRT (HRT), dCRT (dCRT), and our version of CRT (dark green line -CRT-logit).

1 |

 1 1), then step 1 in the procedure defined in Definition B.2 is equivalent with finding k such thatk = max k ∈ [p] | p(k) ≤ k ᾱ . (17)For every i, j, k ∈ [p], let us definep i,j,k = P p i ∈ ((j -1)ᾱ, j ᾱ] , i ∈ S and | S| = k if j ≥ 2 P p i ∈ [0, ᾱ] , i ∈ S and | S| = k if j = 1. (18)Then, since i ∈ S and| S| = k implies that p i ≤ p k ≤ k ᾱ = k ᾱ,we have | S ∩ S c | | S| ∨ 1 = S|=k and i∈ S and 0≤ pi≤k ᾱ .Taking an expectation and writing that 1 0≤ pi≤k ᾱ = 1 pi∈[0, ᾱ]

  Taking the limit where n → ∞ and p fixed, we have, using the result in Theorem 3

Figure 6 :

 6 Figure6: Decoding maps of Relational task in semi-realistic HCP900 experiment, using 400 subjects and dimension reduction to 1000 clusters (i.e. one random seed for generating labels y). We omit Holdout Randomization Test (HRT) as the method does not select any brain region. For dlasso, dCRT and CRT-logit, we plot the test-statistics; for KO the sign of selected coefficients, and for CRT the -log 10 of the empirical p-values.

  Using Eq. (10)

	7	βdy,j ← ( β1, β2, . . . , βj-1, βj+1, . . . , βp)
	8	T decorr j	← decorrelated_test_score(j, X, y,	βd X * ,j , βdy,j ) // Using Eq (11)
	10 end
	11 for j / ∈ ŜSCREENING do
	12	pj = 1
	13 end

9 pj ← 2[1 -|Φ T decorr j |]

Table 1 List of detected genes associated with Glioma Cancer from the TCGA dataset

 1 n = 1026, p = 24776. Empty line (-) signifies no detection. Methods listed in the table are the clustering version. Commonly detected genes between methods are put in bold text. Most detected genes are listed in the mutant list database that can be found in the recorded patients

Table 2 Average runtime of benchmarked methods for one simulation (in seconds).

 2 Standard error is reported in parentheses.

	Methods	Simulated (Sec. 4.2) HCP-semi-real (Sec. 4.3)
	Debiased Lasso (Zhang and Zhang, 2014;	61.83 (5.2)	154.27 (8.79)
	van de Geer et al., 2014; Javanmard and
	Montanari, 2014)
	Knockoff Filter (Barber and Candès, 2015;	1.62 (0.02)	8.12 (0.62)
	Candès et al., 2018)
	CRT (500 samplings) (Candès et al., 2018)	2312.91 (38.21)	7069.96 (109.09)
	HRT (5000 samplings) (Tansey et al.,	14.84 (2.01)	52.17 (4.11)
	2022)
	dCRT[screening=True] (Liu et al., 2021)	16.83 (1.89)	65.18 (3.91)
	dCRT[screening=False] (Liu et al., 2021)	370.12 (8.18)	962.40 (20.63)
	CRT

-logit[screening=True] (this work) 14.16 (0.35) 61.26 (3.55) CRT-logit[screening=False] (this work) 367.91 (4.11) 983.78 (17.26)

  and extended further in[START_REF] Ostrovskii | Finite-sample analysis of m-estimators using self-concordance[END_REF]), that is, |g ′′ (t 1 ) -g ′′ (t)| ≤ |t 1 -t|g ′′ (t) for a fixed constant t, and for every t 1 ∈ R such that t 1 converges to t, with t 1 = β, and t = β 0 . By Assumption 3.1-A3 that X i,j is sub-exponential, applying Bernstein inequality leads to To bound B 3 , by direct application of Hoeffding inequality, we have B 3 ≾ (s

	B 2 ≾	s * log p n	.

* ∨ s ′ ) log p n .

Acknowledgements

BN, BT and SA acknowledged the support of the French "Agence Nationale de la Recherche" under the project ANR-17-CE23-0011 (FastBig) and ANR-20-CHIA-0025-01 (KARAIB AI chair). BN was also supported by Chair DSAIDIS of Télecom Paris.

Appendix A: Proofs of theoretical results in Section 3

We first present some technical lemmas that are useful for the proof of the main theorem. From now on, let ≾ and ≿ denote inequalities with a hidden constant factor, i.e. x ≾ y means that with high probability, there exists an absolute constant C > 0 such that x ≤ Cy, and vice versa. As mentioned in the main text, in what follows, without writing it explicitly, we consider p = p(n).

Lemma A.1 (Lemma E.1, [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF]). Assume Assumption ??, under the logistic model, we have

where s * = ∥β 0 ∥ 0 . In addition, we also have

where g(x) = 1/(1 + exp(x)) is the sigmoid function.

Lemma A.2 (Lemma E.2 [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF], concentration of the gradient and Hessian of the logistic loss function). Assume Assumption ?? holds, under logistic model, we have, with [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF]). Assume Assumption ?? holds, under logistic model, we have

where s * = ∥β 0 ∥ 0 and s ′ = ∥w 0,j ∥ 0 . In addition, we also have

Lemma A.4 (Lemma E.4, [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF], local smoothness conditions on the loss function). Let βnull = (0, β-j ) ∈ R p , where β is an estimator of β 0 . It holds that

for both β = βnull and β = β.

Proof of Theorem 3.1. The following proof is an adaptation from [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF].

Notice that our version of the proof is shorter, with specific consideration on sparse logistic regression, and with elaboration on the convergence rate of the decorreleated test score, which is missing from [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF].

Denote v def.

= (1, ( βd X * ,j ) ⊤ ), then the decorrelated test score can be written in more general from as

Moreover, denote βnull def. = (0, β-j ) and v * def. = (1, w 0,j ), then we have, under the null hypothesis,

A core issue is the dependency of the statistical power and FDR of CRT-logit on the ℓ 1regularization parameter λ dx when doing Lasso distillation on x j in Eq. ( 10). One might choose the heuristic value λ univ = n -1 log p with theoretical validity, as suggested in [START_REF] Ning | A general theory of hypothesis tests and confidence regions for sparse high dimensional models[END_REF]; van de Geer et al. ( 2014). However, experimental results in Fig. 5 show that at λ dx = λ univ (or log 10 λ/λ univ = 0.0 with the labeling of the figure), we do not have the best possible FDR/Power with CRT-logit inference. For this experiment, we average the inference results of 100 simulations (with similar setting in Section 4.1) for different values of n and λ dx , with p fixed. There is a clear phase transition in both FDR and average power when the regularization parameter λ dx increases. In other words, we have found empirically that both FDR and power of the method are sensitive to the ℓ 1 -regularization parameter. Preferably, one wants to return a high statistical power while controlling FDR under predefined level. Hence, it is necessary to choose λ dx wisely. In a more practical scenario, we advise to use cross-validation for X * ,j -distillation operator, as defined by Eq. ( 10). This means we would have to find p different values of λ dx with cross-validation, and we reemphasize the importance of the screening step to reduce the number of computations.

-2.00 -1.75 -1.50 -1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 log 10 λ/λ univ 1200 800 600 400 200 100 No. Samples (n) 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.011 0.035 0.055 0.072 0.095 0.000 0.000 0.000 0.000 0.000 0.000 0.003 0.017 0.042 0.057 0.076 0.098 0.000 0.000 0.000 0.000 0.000 0.000 0.005 0.026 0.049 0.064 0.082 0.102 0.000 0.000 0.000 0.000 0.000 0.001 0.010 0.034 0.055 0.071 0.093 0.102 0.000 0.000 0.000 0.000 0.000 0.005 0.023 0.049 0.063 0.082 0.101 0.102 0.000 0.000 0.000 0.000 0.001 0.013 0.040 0.055 0.074 0.096 0.102 0.102 FDR, α = 0.1 -2.00 -1.75 -1.50 -1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 log 10 λ/λ univ 0.000 0.000 0.001 0.007 0.023 0.072 0.204 0.413 0.595 0.708 0.753 0.791 0.000 0.000 0.002 0.012 0.033 0.106 0.269 0.486 0.648 0.726 0.767 0.798 0.000 0.001 0.004 0.016 0.044 0.134 0.327 0.532 0.674 0.737 0.780 0.800 0.000 0.001 0.007 0.023 0.070 0.195 0.403 0.587 0.706 0.752 0.790 0.799 0.000 0.002 0.015 0.043 0.128 0.315 0.523 0.670 0.735 0.778 0.800 0.799 0.001 0.008 0.028 0.085 0.232 0.450 0.620 0.717 0.757 0.795 0.799 0.799

Avg. Power

Figure 5: FDR/Average Power of 100 runs of simulations while varying the number of samples and ℓ 1 regularization parameter and fixing the number of variables. Note: λ dx is scaled with the factor λ univ = log(p)/n, e.g. the first value for regularization grid is λ dx = 10 -2 λ univ . Default parameter (similar settings in Section 4.1): p = 400, SNR=3.0 (signal-to-noise ration), ρ = 0.5 (feature correlation), κ = 0.05 (sparsity). FDR is controlled at level α = 0.1.

Appendix D: Pseudocode for CRT-logit and Related Algorithms Appendix E: Time complexity of Related Methods

We present the time complexity of benchmarked methods in Table 3.

Algorithm 2: Conditional Randomization Test [START_REF] Candès | Panning for gold: 'model-x' knockoffs for high dimensional controlled variable selection[END_REF] 1 INPUT dataset (X, y), with X ∈ R n×p , y ∈ R n , number of sampling runs B, test statistic Tj, conditional distribution P j|-j for each j = 1, . . . , p ; 2 OUTPUT vector of p-values {pj} p j=1 ; 3 for j = 1, 2, . . . , p do 4 Compute test statistics Tj for original variable; • βdy,j (λ) ← solve_sparse_logistic_cv(X-j, y) // Using Eq. (2)

• βd X * ,j (λ) = argmin β∈R p-1 1 2 ∥X * ,j -X-jβ∥ 2 2 + λ dx ∥β∥ 1 // with λ dx set using cross-validation 2. Obtain test statistic:

Appendix F: Additional Details on Experiments in Section 4

F.1. Preprocessing of the brain-imaging dataset

The Human Connectome Project dataset (HCP) is a collection of brain imaging data on healthy young adult subjects with age ranging from 22 to 35. The participants performed different tasks while being scanned by a magnetic resonance imaging (MRI) device to record blood oxygenation level dependent (BOLD) signals of the brain. The aim of this analysis is to investigate which areas of the brain can predict cognitive activity across participants, while taking into account the information from other brain regions. The brain imaging modalities include, among others, resting-state fMRI (R-fMRI) and taskevoked fMRI (T-fMRI). In this work, we only deal with decoding the task-evoked fMRI dataset. The four classification problems we are working with are as follows.

• Relational: predict whether the participant matches figures or identified feature similarities. • Gambling: predict whether the participant gains or loses gambles.

Algorithm 4: Holdout Randomization Test [START_REF] Tansey | The Holdout Randomization Test for Feature Selection in Black Box Models[END_REF] 1 INPUT dataset (X, y), with X ∈ R n×p , y ∈ R n , number of sampling runs B, test statistic Tj, conditional distribution P j|-j for each j = 1, . . . , p, empirical risk L(•) ; 2 OUTPUT vector of p-values {pj} p j=1 ; 3 (Xtrain, ytrain), (Xtest, ytest) ← data_splitting(X, y); 4 fθ ← model_fitting(Xtrain, ytrain); 5 for j = 1, 2, . . . , p do 6 Tj ← L(Xtest, ytest, fθ (Xtest)); • Emotion: predict whether the participant watches an angry face or a geometric shape. • Social: predict whether the participant watches a movie with social behavior or not.

To perform dimension reduction, our goal is to apply a clustering scheme that keeps the spatial structure of the data. This is achieved with data-driven parcellation along with a spatially constrained clustering algorithm, following the conclusions by [START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF] and [START_REF] Thirion | Which fMRI clustering gives good brain parcellations[END_REF]. The hierarchical clustering scheme that we use recursively merges pair of clusters of features based on a criterion that minimized the within-cluster variance. This algorithm is implemented in scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF], a popular package for applied machine learning.

Appendix G: Example of decoding maps in semi-realistic brain-analysis experiment of Section 4.3