Zipper: Addressing Degeneracy in Algorithm-Agnostic Inference

Authors: geng  chen
yinxu  jia
guanghui  wang
changliang  zou
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 05:33:49

Abstract: The widespread use of black box prediction methods has sparked an increasing interest in algorithm/model-agnostic approaches for quantifying goodness-of-fit, with direct ties to specification testing, model selection and variable importance assessment. A commonly used framework involves defining a predictiveness criterion, applying a cross-fitting procedure to estimate the predictiveness, and utilizing the difference in estimated predictiveness between two models as the test statistic. However, even after standardization, the test statistic typically fails to converge to a non-degenerate distribution under the null hypothesis of equal goodness, leading to what is known as the degeneracy issue. To addresses this degeneracy issue, we present a simple yet effective device, Zipper. It draws inspiration from the strategy of additional splitting of testing data, but encourages an overlap between two testing data splits in predictiveness evaluation. Zipper binds together the two overlapping splits using a slider parameter that controls the proportion of overlap. Our proposed test statistic follows an asymptotically normal distribution under the null hypothesis for any fixed slider value, guaranteeing valid size control while enhancing power by effective data reuse. Finite-sample experiments demonstrate that our procedure, with a simple choice of the slider, works well across a wide range of settings.

Main paper: 
 introduction: 
 Consider predicting response Y ∈ R from covariates X ∈ R p . Due to the popularity of black box prediction methods like random forests and deep neural networks, there has been a growing interest in the so-called "algorithm (or model)-agnostic" inference on the goodness-of-fit (GoF) in regression [1, 2, 3, 4, 5] . This framework aims to assess the appropriateness of a given model for prediction compared to a potentially more complex (often higher-dimensional) model. A common approach involves defining a predictiveness criterion, employing either the sample-splitting or crossfitting strategy to estimate the predictiveness of the two models, and examining the difference in predictiveness. The main focus of this work is to address the issue of degeneracy encountered in predictiveness-comparison-based test statistics.Let P represent the joint distribution of (Y, X) and consider a class F of prediction functions that effectively map the covariates to the response. Define a criterion C( f , P ), which quantifies the predictive capability of a prediction function f ∈ F. Larger values of C( f , P ) indicate stronger predictive capability. The optimal prediction function within the class F is determined as f ∈ arg max f ∈F C( f , P ). Therefore, C(f, P ) represents the highest achievable level of predictiveness within the class F. Examples of C include the (negative) squared loss C( f , P ) = -E[{Y -f (X)} 2 ] for continuous responses and the (negative) cross-entropy loss C( f , P ) = E[Y log f (X) + (1 -Y ) log{1 -f (X)}] for binary responses, where E denotes the expectation under P .In GoF testing problems, there are typically two classes of prediction functions F and F S , where F S is a subset of F. Define a dissimilarity measure ψ S = C(f, P ) -C(f S , P ), which quantifies the deterioration in predictive capability resulting from constraining the model class to F S . Here, f S ∈ arg max f ∈F S C( f , P ) denotes the optimal prediction function within the restricted class. Since F S ⊆ F, it follows that ψ S ≥ 0. The GoF testing is then formulated as H 0 : ψ S = 0 versus H 1 : ψ S > 0.(EQUATIONThe formulation of assessing a scalar predictiveness quantity is inspired by the work of Williamson et al. [5] , which focuses on evaluating variable importance. Other predictiveness or risk quantities have also been explored by [6, 7, 8, 9] . Here, we highlight a few examples where the aforementioned framework can be directly applied.• (Specification testing) Model specification testing aims to evaluate the adequacy of a class of postulated models, such as parametric models, say, examining whether the conditional expectation E(Y | X) = g θ (X) holds almost surely [10] , where g θ is a known function up to an unknown parameter θ. Under the framework (1), we can consider F as a generally unrestricted class, while F S represents a class of parametric models. • (Model selection) GoF testing can also be employed to identify the superior predictive model from a set of candidates [11, 12] . This situation often arises when choosing between two prediction strategies, such as an unregularized model and a regularized one. Testing H 0 is to assess whether the inclusion of a regularizer provides benefits for predictions. • (Variable importance measure) A recently popular aspect of GoF testing involves evaluating the significance of a specific group of covariates U in predicting the response, where X = (U ⊤ , V ⊤ ) ⊤ . This assessment can be incorporated into (1) by defining a subset of prediction functions F S ⊆ F that disregard the covariate group U when making predictions. When utilizing the squared loss, ψ S simplifies to the LOCO (Leave Out COvariates) variable importance measure [1, 13, 14, 5, 15, 16] , which quantifies the increase in prediction error resulting from the removal of U . Specifically, takingψ S = E[{Y -E(Y | V )} 2 ] -E[{Y - E(Y | X)} 2] corresponds to testing for conditional mean independence [17, 18] .The measure ψ S possesses the advantages of being both model-free and algorithm-agnostic. It is not tied to a specific model and remains independent of any particular model fitting algorithm. This flexibility makes ψ S a versatile quantity for evaluating goodness-of-fit, enabling us to leverage diverse machine learning prediction algorithms in its estimation.Let Z = (Y, X) ∼ P , and suppose we have collected a set of independent realizations of Z as Z i = (Y i , X i ) for i = 1, . . . , n. To effectively estimate ψ S while ensuring algorithm-agnostic inference, the sample-splitting or cross-fitting has recently gained significant popularity. This approach relaxes the requirements imposed on estimation algorithms, allowing for the utilization of flexible machine learning techniques [13, 19, 2, 20, 21] . Taking sample-splitting as an example, the data is divided into a training set and a testing set. Based on the training data, estimators f tr n and f tr n,S are obtained for the optimal prediction functions f and f S , respectively. The dissimilarity measure ψ S can then be evaluated using the testing data, i.e., ψ n,S = C(f tr n , P te n ) -C(f tr n,S , P te n ), where P te n represents the empirical distribution function of the testing data. Notably, in the context of measuring variable importance, it has been established that when ψ S > 0, the estimator ψ n,S exhibits asymptotic linearity under certain assumptions [5] . Similar results are also applicable to the problem of comparing multiple algorithms or models [11, 12] . However, the situation becomes distinct when considering the null hypothesis H 0 : ψ S = 0. Recent studies have drawn significant attention to a phenomenon known as degeneracy [22, 23, 15, 18, 5] .Consider the simplest scenario where there are no covariates and the objective is to test whether µ := E(Y ) = 0; see also Example 1 in Lei [11] . In this case, we set F = R and F S = {0}. Using the squared loss, we haveψ S = E(Y 2 ) -E{(Y -µ) 2 } = µ 2 . The estimator for ψ S based on sample-splitting is ψ n,S = 2 Ȳ te n Ȳ tr n -( Ȳ tr n ) 2, where Ȳ tr n and Ȳ te n denote the sample means of the training and testing data, respectively. When µ ̸ = 0, √ n(ψ n,Sµ) follows an asymptotic normal distribution. However, when µ = 0, √ nψ n,S = O P (n -1/2 ), indicating the presence of the degeneracy phenomenon. While inference at a n-rate remains feasible under degeneracy in this specific example, it is crucial to recognize that this would not hold true for intricate models and black box fitting algorithms. Williamson et al. [5] provide evidence for the occurrence of degeneracy in a general variable importance measure ψ S , where the influence function becomes exactly zero under the null hypothesis. It is therefore required to address this degeneracy issue in a generic manner to perform algorithm-agnostic statistical inference.By using the fact that the influence functions of the individual components C(f, P ) and C(f S , P ) in ψ S remain nondegenerate even under H 0 , Williamson et al. [5] proposed an additional data split of the testing data, where the two predictiveness functions are evaluated on two separate testing data splits. This approach ensures a nondegenerate influence function under H 0 , therefore restoring asymptotic normality. A similar approach has been independently explored by Dai et al. [23] . However, performing additional data splits may reduce the actual sample size used in the testing, leading to a substantial loss of power. Alternatively, Rinaldo et al. [13] and Dai et al. [23] introduced data perturbation methods, where independent zero-mean noises are added to the empirical influence functions. These methods also restore asymptotic normality. However, determining the appropriate amount of perturbation to achieve desirable Type I error control remains a heuristic process. Furthermore, Verdinelli and Wasserman [22] proposed expanding the standard error of the estimator to mitigate the impact of degeneracy.In this paper, we introduce the Zipper device for algorithm-agnostic inference under the null hypothesis H 0 of equal goodness. Our approach is inspired by the method of additional splitting of testing data, as demonstrated in the works of Williamson et al. [5] and Dai et al. [23] for assessing variables with zero-importance. Instead of creating two distinct testing data splits to evaluate the discrepancy of predictiveness criteria between the expansive and restricted models, we encourage an overlap between them. The Zipper device serves to bind together the two overlapping splits, with a slider parameter τ ∈ [0, 1) controlling the proportion of overlap. To ensure stable inference and accommodate versatile machine learning prediction algorithms, we incorporate a K-fold cross-fitting scheme.We will demonstrate that the proposed test statistic follows an asymptotically normal distribution H 0 for any fixed value of τ , ensuring valid size control while providing satisfactory power enhancement.For variable importance assessments, in addition to LOCO methods within our framework, Shapley value-based measures are commonly used [24, 25, 26, 27] . These measures, which estimate the incremental predictive accuracy contributed by a specific variable across all possible covariate subsets, reveal complex inter-variable relationships but at a considerable computational expense. Furthermore, conditional randomization tests [28, 29] offer a robust alternative when covariate distributions are known or can be accurately estimated. These methods are especially beneficial in semi-supervised settings with extensive unlabeled data. Additionally, LIME [30] focuses on estimating variable importance locally for a specific instance, while LOCO methods are designed to assess global variable importance.The remainder of our paper is structured as follows. In Section 2, we introduce the Zipper device for addressing the degenerate issue, and present the asymptotic behaviors of the method. Finite-sample experiments are presented in Section 3. Section 4 concludes the paper. Proofs of theorems and additional numerical results are provided in Appendix.2 Our Remedy 
the zipper device: 
 To initiate the process, we randomly partition the data into K folds, denoted as D 1 , . . . , D K , ensuring that each fold is approximately of equal size. For a given fold index k ∈ {1, . . . , K}, we construct estimators f k,n and f k,n,S for the oracle prediction functions f and f S correspondingly, using the data excluding the fold D k . To activate the Zipper device, we further randomly divide D k into two splits, labeled as D k,A and D k,B , with approximately equal sizes, allowing for an overlap denoted as D k,o . Specifically, D k,A ∪ D k,B = D k and D k,A ∩ D k,B = D k,o . Let D k,a = D k,A \D k,C k,n = τ C(f k,n , P k,n,o ) + (1 -τ )C(f k,n , P k,n,a ) and C k,n,S = τ C(f k,n,S , P k,n,o ) + (1 -τ )C(f k,n,S , P k,n,b )represent weighted aggregations of empirical predictiveness criteria corresponding to overlapping and non-overlapping parts, used for estimating C(f, P ) and C(f S , P ), respectively. By employing the cross-fitting process, we obtain the final estimator of ψ S , denoted asEQUATIONwhich is the average over all testing folds. Notably, if we fully open the Zipper, setting τ = 0 and leaving D k,o empty, our approach aligns with the vanilla data splitting strategy utilized in Williamson et al. [5] and Dai et al. [23] for assessing variables with zeroimportance. Conversely, when we completely close Zipper with τ = 1, our method essentially involves evaluating the predictiveness discrepancy using the identical data splitD k,o = D k,A = D k,B = D k ,which is known to result in the phenomenon of degeneracy [5] . Therefore, to avoid such degeneracy, we restrict the slider parameter τ ∈ [0, 1); see also Remark 2.5 below.We start by investigating the asymptotic linearity of the proposed test statistic ψ n,S in (2), which serves as a foundation for establishing the asymptotic distribution under the null hypothesis, as well as for analyzing the test's power. Theorem 2.1 (Asymptotic linearity). If Conditions (C1)-(C4) in Section A hold for both tuples (P, F, f, f k,n ) and (P, F S , f S , f k,n,S ), thenEQUATIONwhere ϕ(Z) = Ċ(f, P ; δ Z -P ) and ϕ S (Z) = Ċ(f S , P ; δ Z -P ). Here, Ċ( f , P ; h) represents the Gâteaux derivative of P → C( f , P ) at P in the direction h, and δ z denotes the Dirac measure at z. Consequently, for any τ ∈ [0, 1), [5] . Given the asymptotic linearity of ψ n,S , we can readily obtain its asymptotic distribution by observing the independence of data in{n/(2 -τ )} 1/2 (ψ n,S -ψ S ) d → N (0, ν 2 S,τ ) as n → ∞, where ν 2 S,τ = (1 -τ )(σ 2 + σ 2 S ) + τ η 2 S , σ 2 = E{ϕ 2 (Z)}, σ 2 S =∪ K k=1 D k,a , ∪ K k=1 D k,b and ∪ K k=1 D k,o , as well as the facts that | ∪ K k=1 D k,a | = | ∪ K k=1 D k,b | = (1 -τ )n/(2 -τ ) and | ∪ K k=1 D k,o | = τ n/(2 -τ ).The asymptotic linearity of ψ n,S in (3) exhibits distinct behaviors depending on the validity of H 0 . Under H 0 : ψ S = 0, we have ϕ = ϕ S almost surely, due to f = f S almost surely. Consequently, the overlapping terms i:Zi∈D k,o {ϕ(Z i )ϕ S (Z i )} for k = 1, . . . , K vanish. As a result, (3) simplifies toψ n,S -ψ S = 1 n/(2 -τ ) K k=1 i:Zi∈D k,a ϕ(Z i ) - K k=1 i:Zi∈D k,b ϕ S (Z i ) + o P (n -1/2 ). (4)In this case, the dominant term compares the means of two independent samples, each with a size of (1τ )n/(2τ ). Additionally, it can be deduced that ν 2 S,τ = (1τ )(σ 2 + σ 2 S ) under H 0 . Conversely, under H 1 : ψ S > 0, we can reformulate (3) asψ n,S -ψ S = 1 n/(2 -τ ) K k=1 i:Zi∈D k,A ϕ(Z i ) - K k=1 i:Zi∈D k,B ϕ S (Z i ) + o P (n -1/2 ).Here, the leading term compares the means of two samples that have an overlap, with each sample having a size of n/(2τ ).To conduct the test, it is crucial to have a consistent estimator for the variance ν 2 S,τ = (1τ )(σ 2 + σ 2 S ) under H 0 . Following the plug-in principle, we derive an estimator ν 2 n,S,τ := (1 -τ )K -1 K k=1 (σ 2 k,n + σ 2 k,n,S ), where for each k ∈ {1, . . . , K}, In many cases, the predictiveness measure takes the form C( f , P ) = E P {g(Y, f (X))}, where g is a known function. In such situations, it has been found that the Gâteaux derivative can be expressed as ϕ f (z) := g(y, f (x)) -E{g(Y, f (X))}; see, for example, Appendix A in Williamson et al. [5] . Consequently, when f = f , the empirical Gâteaux derivative becomes Ċ(f k,n , P k,n ; σ 2 k,n = 1 |D k | i:Zi∈D k { Ċ(f k,n , P k,n ; δ Zi -P k,n )} 2 , σ 2 k,n,S = 1 |D k | i:Zi∈D k { Ċ(f k,n,S , P k,n ; δ Zi -P k,n )}δ z -P k,n ) = g(y, f k,n (x)) -n -1 k i:Zi∈D k g(Y i , f k,n (X i )).T τ := {n/(2 -τ )} 1/2 ψ n,S /ν n,S,τ d → N (0, 1)under H 0 for any τ ∈ [0, 1). For a prespecified significance level α ∈ (0, 1), we reject the null hypothesis if T τ > z 1-α , where z α denotes the α quantile of N (0, 1). A summary of the entire testing procedure can be found in Section C. Remark 2.5. In our asymptotic analysis of the null behavior, we explicitly exclude the case of τ = 1 due to the resulting degeneracy phenomenon. Specifically, under H 0 , when τ = 1, the linear leading term of (3) becomes exactly zero. Therefore, including τ = 1 can introduce a distortion in the Type I error.Next, we turn our attention to the power analysis of the proposed test under the alternative hypothesisH 1 : ψ S > 0.Theorem 2.6 (Power approximation). Suppose the conditions stated in Theorem 2.1 and Proposition 2.3 hold. Then for any τ ∈ [0, 1), the power functionEQUATIONS,τ = {(1τ )(σ 2 + σ 2 S )} 1/2 and Φ denotes the distribution function of N (0, 1). Furthermore, if Cov{ϕ(Z), ϕ S (Z)} ≥ 0, then G S,n,α (τ ) increases with τ .Remark 2.7. The form of the power function can be directly derived from Theorem 2.1 and the fact that the estimator of standard deviation ν n,S,τ p → ν(0) S,τ as n → ∞ under H 1 . Recall that ν 2 S,τ = (1 -τ )(σ 2 + σ 2 S ) + τ η 2 S, which is provably a decreasing function of τ when Cov{ϕ(Z), ϕ S (Z)} ≥ 0. Consequently, the approximate power function G S,n,α (τ ) increase with τ . For more details, please refer to Section B.3. Remark 2.8. The requirement Cov{ϕ(Z), ϕ S (Z)} ≥ 0 is relatively benign. For instance, when consideringψ S = E[{Y -E(Y | V )} 2 ] -E[{Y -E(Y | X)} 2 ]within the framework of evaluating variable importance, this condition is readily satisfied; refer to B.3.1.Consider the "unzipped" version of the proposed test with τ = 0, which has been explored in the works of Williamson et al. [5] and Dai et al. [23] for assessing variable importance. According to Theorem 2.6, the approximate power function reduces toG S,n,α (0) = Φ -z 1-α + (n/2) 1/2 ψ S (σ 2 + σ 2 S ) 1/2 ,aligning with the findings in Dai et al. [23] . In contrast, for the Zipper with τ ∈ [0, 1), we haveG S,n,α (τ ) (i) ≥ Φ -z 1-α + {n/(2 -τ )} 1/2 ψ S (σ 2 + σ 2 S ) 1/2 (ii) ≥ G S,n,α (0),due to the facts that ν(0) S,τ ≤ ν S,τ , and ν S,τ ≤ (σ 2 +σ 2 S ) 1/2 if Cov{ϕ(Z), ϕ S (Z)} ≥ 0.Consequently, our method surpasses the vanilla sample-splitting or cross-fitting based inferential procedures that correspond to τ = 0. The improved power can be attributed to two sources: the introduction of the overlap mechanism τ (corresponding to Inequality (ii)), and the utilization of the variance estimator ν 2 n,S,τ (Inequality (i)).Remark 2.9. As discussed, ν 2 n,S,τ is inconsistent for the limiting variance ν 2 S,τ under H 1 when 0 < τ < 1. If the objective is to construct a valid confidence interval for the dissimilarity measure ψ S , it is crucial to use a consistent variance estimator regardless of whether H 0 holds or not. This can be achieved by incorporating an additional plug-in estimator of η 2 S , which is a component of the asymptotic variance ν 2 S,τ . For detailed construction of a valid confidence interval, please refer to Section D in Appendix.Our asymptotic analysis demonstrates that the Zipper device ensures a valid testing size for any fixed slider parameter τ ∈ [0, 1). Moreover, as the slider τ moves away from 0, the power improves. In practical scenarios with finite sample sizes, selecting an appropriate value of τ involves a tradeoff between efficiency and degeneracy. Opting for a larger value of τ can indeed enhance the testing power. However, an excessively large τ approaching 1 would result in degeneracy and potential size inflation. This occurs because the normal approximation (4) breaks down under the null hypothesis. It is worth emphasizing that using a relatively small τ is generally safer and yields improved power compared to the vanilla splitting-based strategies with τ = 0.To achieve better power while maintaining a reliable size, we propose a simple approach for selecting τ . By (4), the asymptotic normality relies on comparing means from two independent samples∪ K k=1 D k,a and ∪ K k=1 D k,b , each with a size of (1 -τ )n/(2 -τ ).To ensure a favorable normal approximation, we can choose the sample size (1τ )n/(2τ ) such that it meets a predetermined "large" sample size, such as n 0 = 30 or 50. Say, we can specify τ = τ 0 := (n -2n 0 )/(nn 0 ). In the case of very large n, a truncation may be needed to safeguard against degeneracy. For example, we can set τ = min{τ 0 , 0.9}. Our numerical experiments show that this selection of τ achieves satisfactory performances across a wide range of scenarios. For more details on the impact of τ , please refer to Section E.1.3 Finite-Sample Experiments 3.1 Synthetic ExperimentsQuestion: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: We discuss our limitations and further developments in Section 4.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
variable importance assessment: 
 For illustration, we conduct a simulation study to evaluate the performance of the proposed Zipper method in assessing variables with zero-importance, an area of active research. We compare the empirical size and power of Zipper against several benchmark procedures. Firstly, we consider Algorithm 3 proposed in Williamson et al. [5] (referred to as WGSC-3) and the two-split test in Dai et al. [23] (DSP-Split). Both procedures involve an additional splitting of the testing data and can be seen as approximate counterparts to the proposed Zipper test with τ = 0. Another benchmark procedure is Algorithm 2 from Williamson et al. [5] (WGSC-2), which can be viewed as a rough equivalent to Zipper with τ = 1. Additionally, we include the data perturbation method proposed by Dai et al. [23] (DSP-Pert). For each benchmark procedure, we follow the suggestions of the respective authors to select nuisance parameters. We specify the slider parameter τ = min{τ 0 , 0.9} with n 0 = 50 as suggested in Section 2.5.We consider two models: one with a normal response Y ∼ N (X ⊤ β, σ 2 Y ), and another with a binomial response Y ∼ binom(1, logit(X ⊤ β)), where logit(t) = 1/{1 + exp(-t)}. Both models assume that X ∼ N (0, Σ), where Σ = (0.2 |i-j| ) p×p . For each model, we examine two scenarios. The first scenario is a low-dimensional setting with p ∈ {5, 10} and β = (δ, δ, 5, 0, 5, 0 p-5 ) ⊤ . The second scenario is a high-dimensional setting with p ∈ {200, 1000} and β = (δ, δ, 5 0.01p , 0 ⊤ 0.99p-2 ) ⊤ , where a q represents a q-dimensional vector with all entries set to a. In the normal model, we specify σ 2 Y such that the signal-to-noise ratio β ⊤ Σβ/σ 2 Y = 3 by assigning δ = 0. The objective is to test whether the first two variables contribute significantly to predictions from a sequence of n ∈ {200, 500, 1000} independent realizations of (Y, X).Let F represent a generally unrestricted model class, subjecting to a degree of sparsity under the high-dimensional settings. Consider F S ⊆ F such that the prediction functions exclude the first two components of the covariates. To test the irrelevance of the first two variables in predictions, we examine H 0 : ψ S = 0 in (1). We adopt the squared loss for normal responses and the cross-entropy loss for binomial responses. The ordinary least-squares regression and the LASSO are utilized under the low-dimensional and high-dimensional scenarios, respectively. The significance level is chosen as α = 5%, and our experiments entail 1, 000 replications. These experiments are executed on an Intel Xeon Gold 5118 CPU @ 2.30GHz.Table 1 displays the empirical sizes for different testing procedures with n = 500 and p ∈ {5, 1000}.The results reveal that the Zipper, WGSC-3, and DSP-Split consistently maintain the correct size across all models, as anticipated. In contrast, the WGSC-2 exhibits conservative behavior in the low-dimensional setting and inflated sizes in the high-dimensional settings, primarily due to the degeneracy phenomenon. In addition, the data perturbation method, DSP-Pert, fails to control the size in some cases, particularly in the high-dimensional settings. This instability can be attributed to the selection of the amount of perturbation. Figure 2 depicts the empirical power of various testing methods as a function of the magnitude δ representing variable relevance, when n = 500 and p ∈ {5, 1000}. As expected, the Zipper shows a substantial improvement in power compared to the vanilla crossfitting based approaches, WGSC-3 and DSP-Split, with WGSC-3 and DSP-Split demonstrating similar performances. Under the high-dimensional settings, the WGSC-2 and DSP-Pert exhibit higher power than Zipper, but at the expense of losing valid size control. For a comprehensive analysis of the empirical sizes and power across various combinations of n and p, please refer to Section E.2. These additional results consistently support the conclusion that the Zipper method demonstrates reliable empirical size performance and significant power enhancement compared to that methods that utilize non-overlapping splits.We extend our investigation beyond the variable importance assessment problem to include an evaluation of the proposed Zipper device in addressing model specification issues. Our analysis focuses on the model defined as Y = Xβ + ε, where X ∼ N (0, Σ) with Σ = (0.2 |i-j| ) p×p and ε ∼ N (0, σ 2 ε ). Here, we assume that ∥β∥ 0 = 2, indicating the presence of two nonzero components in β, while the positions of these components remain unknown. Our objective is to test the model specification hypothesis: H 0 : β = ( * , * , 0 p-2 ) ⊤ versus H 1 : ∥β∥ 0 = 2 but not H 0 , where * represents any nonzero value. To generate the data, we consider three scenarios: (i) β = (0.4, 0.4, 0 p-2 ) ⊤ , (ii) β = (0.4, 0, 0.4, 0 p-3 ) ⊤ , and (iii) β = (0, 0, 0.4, 0.4, 0 p-4 ) ⊤ . We determine the value of σ 2 ε such that the signal-to-noise ratio β ⊤ Σβ/σ 2 ε = 1. Subsequently, we generate independent realizations (Y i , X i ) for i = 1, . . . , n, with n = 500 and p ∈ {5, 1000}. To estimate β, we employ ordinary least squares under H 0 , and perform the best two subset selection under H 1 . For the case when p = 5, we conduct an exhaust search. For the case when p = 1000, we utilize the abess algorithm [31] to approximate the solutions. The results are summarized in Table 2 , where we observe that under Scenario (i), corresponding to H 0 , Zipper maintains correct size control. Furthermore, under Scenarios (ii) and (iii), corresponding to H 1 , Zipper exhibits substantial improvements in power when compared to the vanilla cross-fitting based approaches, WGSC-3 and DSP-Split. We apply the Zipper method to the widely used MNIST handwritten digit dataset [32] . The MNIST dataset consists of size-normalized and center-aligned handwritten digit images, each represented as a 28 × 28 pixel grid (resulting in p = 28 2 = 784). For our analysis, we specifically extract subsets of the dataset representing the digits 7 and 9, following Dai et al. [23] , resulting in a total of n = 14251 images. In Figure 3 , we calculate and graphically represent the average grayscale pixel values for images sharing the same numerals. We divide each image into nine distinct regions, as shown by the blank squares in Figure 3 , with the objective of detecting regions that can effectively distinguish between these two digits. To achieve this, we perform a sequence of variable importance testing to assess the relevance of each region in making predictions while considering the remaining regions. Given the nature of the data, we employ a Convolutional Neural Network (CNN) as the underlying model, leveraging its proven effectiveness in image analysis. In the Zipper approach, we select the slider parameter τ such that n 0 = 50, as recommended in the manuscript. As a benchmark, we adopt WGSC-3 [5] (equivalent to DSP-Split [23] ), which produces valid size, aligning with our approach. We set the predefined significance level for each test as α = 0.05/9, applying the Bonferroni correction to account for multiple comparisons. The discovered regions, highlighted in red, are presented in Figure 3 . Our findings indicate that the Zipper method outperforms WGSC-3 in identifying critical regions with greater efficacy.We expand the application of our Zipper method to the bodyfat dataset [33] , which provides an estimate of body fat percentages obtained through underwater weighing, along with various body circumference measurements from a sample of n = 252 men. Our objective is to conduct marginal variable importance tests for each body circumference while considering potential influences from essential attributes such as age, weight, and height. To accurately estimate the relevant regression functions within this dataset, we employ the random forest as our modeling technique. Table 3 presents the resulting p-values obtained from the Zipper and WGSC-3 methods. By applying the Bonferroni correction, the Zipper method identifies both Abdomen and Hip as significant factors at the significance level of α = 0.05/10. In contrast, WGSC-3 suggests only Abdomen as important. It is worth noting that a recent study by Zhu et al. [34] proposed the formula (Waist+Hip)/Height as a straightforward evaluation index for body fat. Remarkably, our finds align with this fact, further supporting the validity and relevance of our Zipper method in identifying key factors. Assume P ∈ M for a rich class M of distributions. Define the linear space H = {c( P1 -P2 ) :c ≥ 0, P1 , P2 ∈ M}. For each h ∈ H, let ∥h∥ ∞ = sup z | F1 (z) -F2 (z)|,where Fj denotes the distribution function corresponding to Pj for j = 1, 2. Consider F as the class of predictions functions endowed with a norm ∥ • ∥ F . Let f ∈ arg max f ∈F C( f , P ), and {f k,n } K k=1 be a sequence of estimators of f . For each 1 ≤ k ≤ K, define a k,n : z → Ċ(f k,n , P ; δ z -P ) -Ċ(f, P ; δ z -P ) and b k,n : z → Ċ(f k,n , P k,n ; δ z -P k,n )-Ċ(f k,n , P ; δ z -P ). The following conditions are imposed on the tuple (P, F, f, f k,n ).(C1) There exists some constant C > 0 such that, for each sequence f1 , f2 , . . .∈ F such that ∥ fj -f ∥ F → 0, |C( fj , P ) -C(f, P )| ≤ C∥ fj -f ∥ 2F for each j large enough; (C2) There exists some constant δ > 0 such that for each sequence ϵ 1 , ϵ 2 , . . . ∈ R and h, h 1 , h 2 , . . . ∈ H satisfying that ϵ j → 0 and ∥h j -h∥ ∞ → 0, it holds that sup f ∈F :∥ f -f ∥ F <δ C( f , P + ϵ j h j ) -C( f , P ) ϵ j -Ċ( f , P ; h j ) → 0; (C3) For each 1 ≤ k ≤ K, ∥f k,n -f ∥ F = o P (n -1/4 ); (C4) For each 1 ≤ k ≤ K, a 2 k,n dP = o P (1); (C5) For each 1 ≤ k ≤ K, b 2 k,n dP = o P (1).C(f k,n , P k,n,I ) -C(f, P ) = 1 |D k,I | i:Zi∈D k,I ϕ(Z i ) + o P (|D k,I | -1/2 ) and C(f k,n,S , P k,n,I ) -C(f S , P ) = 1 |D k,I | i:Zi∈D k,I ϕ S (Z i ) + o P (|D k,I | -1/2 ).Proof. See Theorem 2 in Williamson et al. [5] .Denote m k = |D k,A | = |D k,B |. Since |D k,a | = |D k,b | = (1 -τ )m k and |D k,o | = τ m k , we have m k = n/{(2 -τ )K} due to |D k,a | + |D k,b | + |D k,o | = n/K. By Lemma B.1, we have C k,n -C(f, P ) = τ {C(f k,n , P k,n,o ) -C(f, P )} + (1 -τ ){C(f k,n , P k,n,a ) -C(f, P )} = 1 m k i:Zi∈D k,A ϕ(Z i ) + o P (m -1/2 k).Since the number of folds K is fixed, we haveEQUATIONSimilarly, we conclude thatEQUATIONCombining ( 5) and ( 6), we obtainψ n,S -ψ S = 1 n/(2 -τ ) K k=1 i:Zi∈D k,a ϕ(Z i ) - i:Zi∈D k,b ϕ S (Z i ) + i:Zi∈D k,o {ϕ(Z i ) -ϕ S (Z i )} + o P (n -1/2 ).By applying the standard central limit theorem,EQUATION{(1 -τ )n/(2 -τ )} 1/2 K K k=1 1 (1 -τ )m k i:Zi∈D k,b ϕ S (Z i ) d → N (0, σ 2 S ),as n → ∞. If ψ S > 0, we have{τ n/(2 -τ )} 1/2 K K k=1 1 τ m k i:Zi∈D k,o {ϕ(Z i ) -ϕ S (Z i )} d → N (0, η 2 S ).Hence,{n/(2 -τ )} 1/2 (ψ n,S -ψ S ) d → N (0, ν 2 S,τ ), where ν 2 S,τ = (1 -τ )(σ 2 + σ 2 S ) + τ η 2 S .It suffices to show that σ 2k,n p → σ 2 as n → ∞. Since Ċ(f k,n , P k,n ; δ z -P k,n ) = Ċ(f k,n , P k,n ; δ z -P k,n ) -Ċ(f k,n , P ; δ z -P )+ Ċ(f k,n , P ; δ z -P ) -Ċ(f, P ; δ z -P ) + Ċ(f, P ; δ z -P ),we haveEQUATIONFor any ϵ > 0, by Markov's inequality and Condition (C5),Pr 1 n k i:Zi∈D k { Ċ(f k,n , P k,n ; δ Zi -P k,n ) -Ċ(f k,n , P ; δ Zi -P )} 2 ≥ ϵ | D -k ≤ 1 ϵ E[{ Ċ(f k,n , P k,n ; δ Zi -P k,n ) -Ċ(f k,n , P ; δ Zi -P )} 2 | D -k ] = 1 ϵ {b k,n (z)} 2 dP (z) p → 0.Then, by the law of total expectation, the first term of ( 7) is o P (1). Similarly, the second term is o P (1). Specifically, by Condition (C4), for any ϵ > 0,Pr 1 n k i:Zi∈D k { Ċ(f k,n , P ; δ Zi -P ) -Ċ(f, P ; δ Zi -P )} 2 ≥ ϵ | D -k ≤ 1 ϵ E[{ Ċ(f k,n , P ; δ Z -P ) -Ċ(f, P ; δ Z -P )} 2 | D -k ] = 1 ϵ {a k,n (z)} 2 dP (z) p → 0.For the third term, by the law of large numbers,1 n k i:Zi∈D k Ċ(f, P ; δ Zi -P ) 2 p → E{ Ċ(f, P ; δ Z -P ) 2 } = σ 2 ,as n → ∞. For the fourth term, we have 1 n k i:Zi∈D k { Ċ(f k,n , P k,n ; δ Zi -P k,n ) -Ċ(f k,n , P ; δ Zi -P )} Ċ(f, P ; δ Zi -P ) ≤ 1 n k i:Zi∈D k { Ċ(f k,n , P k,n ; δ Zi -P k,n ) -Ċ(f k,n , P ; δ Zi -P )} 2 1/2 × 1 n k i:Zi∈D k Ċ(f, P ; δ Zi -P )-P k,n ) = g{y, f k,n (x)} - 1 n k i:Zi∈D k g{Y i , f k,n (X i )}, Ċ(f k,n , P ; δ z -P ) = g{y, f k,n (x)} -E P [g{Y, f k,n (X)}].By noting thatE[{ Ċ(f k,n , P k,n ; δ z -P k,n ) -Ċ(f k,n , P ; δ z -P )} 2 | D -k ] = E   1 n k i:Zi∈D k g{Y i , f k,n (X i )} 2 | D -k   -E[g{Y, f k,n (X)} | D -k ] 2 = 1 n k V ar[g{Y, f k,n (X)} | D -k ] p → 0,the conclusion follows.By Theorem 2.1, for any τ ∈ [0, 1),{n/(2 -τ )} 1/2 (ψ n,S -ψ S ) d → N (0, ν 2 S,τ ) as n → ∞, where ν 2 S,τ = (1 -τ )(σ 2 + σ 2 S ) + τ η 2 S , σ 2 = E{ϕ 2 (Z)}, σ 2 S = E{ϕ 2 S (Z)}, and η 2 S = E[{ϕ(Z) -ϕ S (Z)} 2 ]. By examining the proof of Proposition 2.3, we conclude that ν n,S,τ p → ν (0) S,τ under H 1 . Consequently, the power function isPr(T τ > z 1-α | H 1 ) = Pr ν (0) S,τ ν n,S,τ {n/(2 -τ )} 1/2 (ψ n,S -ψ S ) ν S,τ > ν (0) S,τ ν S,τ z 1-α - {n/(2 -τ )} 1/2 ν (0) S,τ ψ S ν n,S,τ ν S,τ | H 1 = Φ - ν (0) S,τ ν S,τ z 1-α + {n/(2 -τ )} 1/2 ψ S ν S,τ + o(1)= G S,n,α (τ ) + o(1).Given that η S = σ 2 + σ 2 S -2Cov{ϕ(Z), ϕ S (Z)}, the asymptotic variance ν 2 S,τ can be expressed asσ 2 + σ 2 S -2τ Cov{ϕ(Z), ϕ S (Z)}. Let φ 1 (τ ) = (1 -τ )(σ 2 + σ 2 S ) σ 2 + σ 2 S -2τ Cov{ϕ(Z), ϕ S (Z)} and φ 2 (τ ) = {n/(2 -τ )} 1/2 ψ S (σ 2 + σ 2 S -2τ Cov{ϕ(Z), ϕ S (Z)}) 1/2 . Since 2Cov{ϕ(Z), ϕ S (Z)} ≤ σ 2 + σ 2S and τ ∈ [0, 1), it follows that φ 1 (τ ) > 0 and φ 2 (τ ) > 0. When Cov{ϕ(Z), ϕ S (Z)} ≥ 0, we observe that φ 1 (τ ) is monotonically decreasing, and φ 2 (τ ) is monotonically increasing with respect to τ . Therefore, the approximate power function G S,n,α (τ ) increase with τ .Here we demonstrate the validity of the condition Cov{ϕ(Z), ϕ S (Z)} ≥ 0 in assessing variable importance when utilizing the square loss. Consider the scenario where Y = E(Y | X) + ε, with the noise ε assumed to be independent of the covariate X. We seek to evaluate the significance of a specific set of covariates U in predicting the response, whereX = (U ⊤ , V ⊤ ) ⊤ . Assuming E(Y | X) ∈ F and E(Y | V ) ∈ F S , the optimal prediction functions within F and F S are f (X) = E(Y | X) and f S (V ) = E(Y | V ), respectively.Referring to Remark 2.4, the covariance Cov{ϕ(Z), ϕ S (Z)} can be expressed asCov{ϕ(Z), ϕ S (Z)} = Cov[{Y -f (X)} 2 , {Y -f S (V )} 2 ] = V ar(ε 2 ) + Cov[ε 2 , {f (X) -f S (V )} 2 ] + 2Cov[ε 2 , ε{f (X) -f S (V )}].It's evident that the first two terms are non-negative. Regarding the third term, since ε is independent of X and Ef (X) | V = f S (V ), we deduceCov[ε 2 , ε{f (X) -f S (V )}] = E[ε 3 {f (X) -f S (V )}] -E(ε 2 )E[ε{f (X) -f S (V )}] = E(ε 3 )E{f (X) -f S (V )} -E(ε 2 )E(ε)E{f (X) -f S (V )} = E(ε 3 )E[E{f (X) -f S (V ) | V }] = 0.Thus, Cov{ϕ(Z), ϕ S (Z)} ≥ 0 follows.Algorithm 1 The algorithm for the proposed Zipper testing procedure Input: Observed data {Z i } n i=1 , number of folds K, and slider parameter τ ∈ [0, 1)Randomly partition {Z i } n i=1 into K disjoint folds D 1 , . . . , D K for k = 1, . . . , K do Using data in {D j } K j=1 \ D k , construct estimators f k,n of f and f k,n,S of f S Using data in D k , construct estimator P k,n of P Randomly divide D k into two splits D k,A and D k,B of roughly equal size m k , with an overlap D k,o such that |D k,o | = ⌊τ m k ⌋. Let D k,a = D k,A \D k,o and D k,b = D k,B \D k,o . Using data in D k,I , construct estimators P k,n,I for I ∈ {o, a, b} Compute C k,n = τ C(f k,n , P k,n,o ) + (1 -τ )C(f k,n , P k,n,a ) and σ 2 k,n = n -1 k i:Zi∈D k { Ċ(f k,n , P k,n ; δ Zi -P k,n )} 2 , where n k = |D k | Compute C k,n,S = τ C(f k,n,S , P k,n,o ) + (1 -τ )C(f k,n,S , P k,n,b ) and σ 2 k,n,S = n -1 k i:Zi∈D k { Ċ(f k,n,S , P k,n ; δ Zi -P k,n )} 2 end for Compute C n = K -1 K k=1 C k,n , C n,S = K -1 K k=1 C k,n,S and estimator ψ n,S = C n -C n,S of ψ S Compute ν 2 n,S,τ = (1 -τ )K -1 K k=1 (σ 2 k,n + σ 2 k,n,S ) Output: P-value 1 -Φ({n/(2 -τ )} 1/2 ψ n,S /ν n,S,τ ), where Φ denotes the distribution function of N (0, 1)To construct a valid confidence interval for the dissimilarity measure ψ S , it is crucial to use a consistent variance estimator irrespective of whether H 0 holds or not. This can be achieved by incorporating an additional plug-in estimator of η 2 S , which is a component of the asymptotic varianceν 2 S,τ . Let ν 2 n,S,τ = K -1 K k=1 {(1 -τ )(σ 2 k,n + σ 2 k,n,S ) + τ η 2 k,n,S }, where for each k ∈ {1, . . . , K}, η 2 k,n,S = 1 n k i:Zi∈D k { Ċ(f k,n , P k,n ; δ Zi -P k,n ) -Ċ(f k,n,S , P k,n ; δ Zi -P k,n )} 2 .The confidence interval for ψ S can be formulated asψ n,S - z α/2 ν n,S,τ {n/(2 -τ )} 1/2 , ψ n,S + z α/2 ν n,S,τ {n/(2 -τ )} 1/2 ,which has an asymptotic coverage of 1α. To assess the influence of the slider parameter τ on the size and power of the proposed Zipper approach, we conduct a small simulation study using the simulation settings outlined in the manuscript. The results of this study are summarized in Tables 4 5 . Notably, we observe that the empirical sizes are consistently controlled when τ falls within an appropriate range. However, selecting excessively large values of τ led to a slight inflation in size, primarily due to the poor normal approximation under the null hypothesis, as discussed in Section 2.5. Furthermore, the empirical power increases as the value of τ increases, which aligns with our expectations.To provide a comprehensive analysis of the empirical sizes and power across various settings, we conducted additional simulations using the simulation settings outlined in the manuscript. We explored different combinations of the sample size n and dimension p. Moreover, we investigated a model with a heavy-tailed response by considering Y = Xβ + σ Y ε/3, ε ∼ t 3 , where t 3 represents the t-distribution with 3 degrees of freedom. The results for empirical sizes and power with δ = 1 are summarized in Tables 6 and 7 . These additional findings consistently support the conclusion that the Zipper method demonstrates reliable empirical size performance and provides significant power enhancement compared to methods that utilize non-overlapping splits.To evaluate variables with zero-importance, Rinaldo et al. [13] and Dai et al. [23] proposed a data perturbation technique, injecting independent zero-mean noises into empirical influence functions. We describe this method using our own terminology. We begin by partitioning the data randomly into K folds, denoted as D 1 , . . . , D K , ensuring equitable sizing across folds. For each fold index k ∈ {1, . . . , K}, we construct estimators f k,n and f k,n,S for the oracle prediction functions f and f S , correspondingly, utilizing data excluding fold D k . This method employs the following test statistic: where ε i ∼ N (0, 1) for i = 1, . . . , n are independent noise, and ρ represents the perturbation parameter governing the extent of perturbation. It can be established that for any ρ > 0,ϕ n,S,pert,ρ = 1 n K k=1 i:Zi∈D k ϕ(Z i ) -ϕ S (Z i ) + ρε i ,n 1/2 (ϕ n,S,pert,ρ -ψ S ) d → N (0, ν 2 S,pert,ρ ), as n → ∞, where ν 2 S,pert,ρ = η 2 S + ρ 2 , and η 2 S = E[{ϕ(Z) -ϕ S (Z)} 2 ]. Following the plug-in principle, we employ the normalized test statistic T pert,ρ := n 1/2 ϕ n,S,pert,ρ /ν n,S,pert,ρ , which converges in distribution to N (0, 1) under H 0 for any ρ > 0.Here ν 2 n,S,pert,ρ = K -1 K k=1 η 2 k,n,S + ρ 2 , where for each k ∈ {1, . . . , K}, η 2 k,n,S = n k -1 i:Zi∈D k { Ċ(f k,n , P k,n ; δ Zi -P k,n ) -Ċ(f k,n,S , P k,n ; δ Zi -P k,n )} 2 .For a prespecified significance level α ∈ (0, 1), this method rejects H 0 if T pert,ρ > z 1-α . Indeed, a direct correspondence emerges between the slider parameter τ within our Zipper and the perturbation parameter ρ employed in the data perturbation technique, thereby facilitating a comparative analysis of the two methodologies. To elucidate this connection, we demonstrate that, under H 1 : ψ S > 0, the power function corresponding to the data perturbation method isPr(T pert,ρ > z 1-α | H 1 ) = Pr n 1/2 (ϕ n,S,pert,ρ -ψ S ) ν S,pert,ρ > z 1-α - n 1/2 ψ S ν S,pert,ρ | H 1 = Φ -z 1-α + n 1/2 ψ S ν S,pert,ρ + o(1) := G S,n,α,pert (ρ) + o(1).First notice that our Zipper method with τ = 0 (i.e., the vanilla sample-splitting) has approximate powerG S,n,α (0) = Φ -z 1-α + (n/2) 1/2 ψ S (σ 2 + σ 2 S ) 1/2 .We establish an upper bound for the perturbation parameter ρ, namely,ρ 2 ≤ 2(σ 2 + σ 2 S ) -η 2 S , where equality yields G S,n,α,pert (ρ) = G S,n,α (0). Should ρ 2 > 2(σ 2 + σ 2 S ) -η 2 S, it becomes evident that G S,n,α,pert (ρ) ≤ G S,n,α (0), suggesting that the data perturbation technique may even exhibit lower power.The idea is to establish a relationship between τ and ρ such that both methods yield similar power, given their valid size control for fixed τ and ρ. Employing a consistent variance estimator as per Remark 2.9, our Zipper method exhibits an approximate power ofG (CI) S,n,α (τ ) = Φ -z 1-α + {n/(2 -τ )} 1/2 ψ S ν S,τ .S,n,α (τ ) = G S,n,α,pert (ρ), we obtainEQUATIONestablishing a one-to-one correspondence. Under this correspondence,EQUATION, where the latter denotes the approximate power of our Zipper method utilizing the proposed variance estimation scheme. Thus, our Zipper method can be regarded as a data-adaptive perturbation strategy that circumvents external randomization while potentially offering power enhancement owing to the variance estimation scheme.Figure 4 illustrates the empirical power comparison between Zipper and the data perturbation method, where the perturbation parameter ρ = ρ(τ ) is determined by (8) , across various values of τ . This assessment is conducted according to the low-dimensional normal response scenario outlined in Section 3, particularly with p = 5 and n = 500. As anticipated, Zipper and the perturbation method demonstrate the same power when τ = 0, while for τ ∈ (0, 1), Zipper consistently exhibits superior power compared to the perturbation method. NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes] Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: We concretely describe our Zipper algorithm in Section C, and our sythetic data generating procedure in Section 3. The datasets utilized in real-data analysis are sourced from publicly available materials on the internet. We think they are sufficient to reproduce the main experimental results. Guidelines:• The answer NA means that paper does not include experiments requiring code. Justification: Our work conducts a general framework of algorithm/model agnostic inference on goodness-of-fit in regression. This is a fundamental research on statistical methodology and has no direct societal impacts as far as we currently understand. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: The paper poses no such risks. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: The datasets used in real-data examples are all open-source datasets, and we properly cite the original paper that produced the dataset. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA]Justification: The paper does not release new assets.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: The paper Addresses degeneracy issues in algorithm/model-agnostic goodness-of-fit inference such as assessing variable importance.
It broadly applicable scenario (model/algorithm agnostic inference).
Also, it does so with minimal intervention to model fitting / testing (just split your data in a certain way, evaluate on carefully defined subsets, and combine).
The authors clearly addressed all questions by the reviewers during the rebuttal period please update the paper with the clarifications discussed.

Reviews: review: paper_summary: paper_summary: This work presents a test of Performance(f, test_data_1) - Performance(f_subset, test_data_2) where f is the best model in a class F and f_subset is the best model in a subset of F. They focus on how precisely to split your samples between test_data_1 and test_data_2 to avoid degeneracies that can arise when asking questions like whether the performance difference is equal to zero.

The main proposed solution is to take an underlying set of test data D_k, split it into two possibly overlapping pieces D_A and D_B, with overlap D_o, and to compute Performance(f, test_data_1) as a combination of performance on D_o and D_A \ D_o, and likewise to compute Performance(f_subset, test_data_2) as a combination of performance on D_o and D_B \ D_o.

For example, Performance(f, test_data_1) = \tau performance(f, D_o) + (1-\tau) performance(f, D_A \ D_o) . Here picking \tau 0<=\tau<1 avoids the degeneracy.  This idea can be seen as executing a particularly careful and favorable property of sample splitting.

They show that the test statistic, under the null hypothesis of no difference, is asymptotically linear. They also derive a consistent estimator for the variance of the statistic under the null.  They also provide selection criteria for the zipper \tau.

Using these tools they analyze the power/size of the test statistic and run simulations in variable importance, model specification, and more.

As future work, they leave open some ideas around which other hyperparameters to select data-adaptively, how the test works on large scale data, and how to control errors.
strength_weakness: strength_weakness: Nothing major at the moment, may revise before discussion period.
questions: questions: Nothing major at the moment, may revise before discussion period.
limitations: limitations: Yes
paper_summary: paper_summary: The authors propose a new test statistic they call Zipper for algorithm-agnostic inference on goodness-of-fit testing. While previous solutions suffer from a degeneracy issue (i.e., fails to converge to a non-degenerate distribution under the null hypothesis of equal goodness), the proposed test statistic does not (i.e., it converges asymptotically to a normal distribution under the null). It also improves power over other proposed solutions which also tackle the degeneracy issue (due to Zipper’s effective reuse of data). Zipper could have applications in specification testing, model selection and variable importance assessments.
strength_weakness: strength_weakness: 1. The paper could benefit from being more self-contained. For instance:
- The theorems reference conditions in the appendix. While I understand why it makes sense to place such details in the appendix, some informal discussion of the conditions in the main paper would be helpful for exposition.
- The paper consistently references WIlliamson et al. without precisely describing what is proposed in the work (beyond the a high-line mention of “Williamson et al. proposed an additional split of the testing data” in lines 97-98). Adding more discussion up front would help with readability / appreciating the paper’s contribution during the analysis.
2. The paper mentions that this test statistic / goodness-of-fit testing broadly is applicable in many settings but does not compare what is proposed with other methods which are also used in these settings (e.g., other measures of variable importance). It would help to add this discussion to the related work. 
3. The numerical experiments / results could be more comprehensive. Specifically:
- The results for the numerical experiments do not make it easy to see that Zipper is preferred over the other methods; for instance, Table 1 shows that Zipper has better size than WSGC-1 and DSP-Pert, but Figure 2 shows that the latter have better power.
- Table 2 is confusing given the caption mentions both size and power but I only see one number for each method. I think based on the text and inferring from the magnitudes that the first row is size while the second and third row are power, but why not include both for all scenarios? 
- Based on Figure 3, it is not obvious to me that the left column (Zipper) is better than the right column (baseline). Could the authors elaborate?
- The text in 3.2.2 makes it more clear what might be preferable about Zipper, but it is not clear based on Table 3 alone. Perhaps some way of directly showcasing the alignment with a previous study would be helpful?
questions: questions: 1. Lines 231-236 suggest that it is not only the zipper mechanism that enables power improvements, but also the use of a variance estimator. Could the authors expand on the latter, e.g. provide more intuition, incorporate it into the overall story in the paper outside of these few lines?
2. When would one prefer this approach to others in, say, variable importance testing (e.g. conditional randomization tests, Shapley values) and why? I understand that the precise goal is slightly different between the options mentioned, but as this work falls under the interpretability / explainability category, it would be useful to situate it explicitly.
3. Could the authors combine the results of Table 1 and Figure 2 to help visualize how Zipper compares holistically to the other methods?
4. Could the authors provide both power and size results in Table 2?
5. Could the authors explain why Figure 3 shows that the proposed method is preferred over the baseline?
limitations: limitations: The applicability of the method for larger scale settings has not yet been tested. Moreover, for an application such as explanability, GoF testing may not necessarily be the most direct given the test merely tests how well a pre-specified subset performs relative to the full set, rather than locating an important subset.

(As a quick aside, I noticed in the checklist that the authors mention discussing limitations in Section 4, but I do not seem to see any mentions of limitations in that section)
paper_summary: paper_summary: This paper proposes a method for quantifying model-agnostic goodness of fit to allow better comparison between two models / model classes etc. In particular, this paper deals with the problem of degeneracy under the null hypothesis of equal goodness. Prior work has addressed this by splitting the test set into distinct subsets, however, this reduces the sample size significantly. To overcome this, Zipper splits the test set into overlapping test sets and uses the proportion of overlap of the splits to better aggregate the test statistic. This way, Zipper is able to more effectively use limited test data.
strength_weakness: strength_weakness: 1. I'm not an expert in this area, but I'm concerned about the novelty and significance of this contribution. I would urge the authors to further clarify the significance of this change in test data splitting (especially for other expert reviewers).
questions: questions: N/A
limitations: limitations: Novelty / Contribution.

Comments: title: On related works and variance underestimation
comment: I have read the rebuttal and in the light of the other reviewer's comments, I change my score to recommend acceptance for this paper.
Thanks for your reply!

- For W2, could the authors provide more discussion one when to consider this class of methods over the others added to the related works discussion?

Response: 

_Shapley value-based methods_: Shapley value-based measures can be conceptualized as a weighted average of LOCO measures within our framework, accounting for the inclusion of a feature (or a feature group) of interest across all possible subsets (Williamson and Feng, 2020; Verdinelli and Wasserman, 2023; 2024). Inference based on Shapley values can utilize sample-splitting methods, as demonstrated by Williamson and Feng (2020). While extending our Zipper framework to accommodate this sample-splitting inference is feasible, it requires the evaluation of all covariate subsets, entailing significant computational resources due to the necessity of fitting numerous models. Moreover, as noted by Verdinelli and Wasserman (2024), many of these submodels may not contribute meaningfully to the overall measure, particularly when each is weighted equally in the definition. We recommend that if there is confidence in the algorithms' ability to approximate true submodels accurately, LOCO measures should be considered, where our framework can be effectively implemented.

_Conditional randomization tests_: CRTs offer flexibility in selecting variable importance measures beyond LOCOs. This method necessitates knowledge of the conditional distribution $X_\mathcal{S}\mid X_{-\mathcal{S}}$, from which data can be sampled for multiple model refitting and measure calculation (Candes et al., 2018). Subsequent research, such as that by Tansey et al. (2021), has aimed to mitigate these computational demands. However, when the distribution of covariates is either unknown or cannot be precisely estimated, our methodology presents a more viable alternative to CRTs.

- For your response Q1, you mentioned, "Under H1, it (the variance estimator) underestimates the actual unknown variance, which is advantageous for power enhancement as demonstrated." Could you explain why that is and where it is demonstrated?

Response: 

_Variance underestimation under $H_1$_: As outlined in Lines 167 and 218, the unknown variance $\nu_{\mathcal{S},\tau}^2=(\nu_{\mathcal{S},\tau}^{(0)})^2+\tau\eta^2_\mathcal{S}$, where $(\nu_{\mathcal{S},\tau}^{(0)})^2=(1-\tau)(\sigma^2+\sigma_{\mathcal{S}}^2)$. Given that $\eta^2_\mathcal{S}\ge 0$ (with $\eta^2_\mathcal{S}=0$ under $H_0$ and $\eta^2_\mathcal{S}>0$ under $H_1$), it follows that $\nu_{\mathcal{S},\tau}^2\ge(\nu_{\mathcal{S},\tau}^{(0)})^2$. Proposition 2.3 and Line 221 confirm that our variance estimator is always consistent to $(\nu_{\mathcal{S},\tau}^{(0)})^2$ under both hypotheses. Thus, under $H_1$, the estimator tends to underestimate $\nu_{\mathcal{S},\tau}^2$.

_Power enhancement_: The aforementioned underestimation leads to an advantageous adjustment in our approximated power function, as outlined in Theorem 2.6 and Line 232. Specifically,

$$
G_{\mathcal{S},n,\alpha}(\tau) = \Phi\left(-\frac{\nu_{\mathcal{S},\tau}^{(0)}}{\nu_{\mathcal{S},\tau}}z_{1-\alpha} + \frac{\{n/(2-\tau)\}^{1/2}\psi_\mathcal{S}}{\nu_{\mathcal{S},\tau}}\right)\ge \Phi\left(-z_{1-\alpha} + \frac{\{n/(2-\tau)\}^{1/2}\psi_\mathcal{S}}{\nu_{\mathcal{S},\tau}}\right),
$$

This expression establishes a lower bound that corresponds to the power function of a test statistic employing a consistent variance estimator for $\nu_{\mathcal{S},\tau}^2$ under both $H_0$ and $H_1$. This explains how the underestimation enhance the test's power. Between Lines 231-236, we further elucidate that this lower bound outperforms the power achievable by conventional sample-splitting-based statistics, thereby underscoring the efficacy of our overlapping mechanism in boosting power.

Thank you once again for your insightful questions.

References

Brian Williamson and Jean Feng. Efficient nonparametric statistical inference on population feature importance using Shapley values. *Proceedings of the 37th International Conference on Machine Learning*, PMLR 119:10282-10291, 2020.

Isabella Verdinelli and Larry Wasserman. Feature importance: A closer look at shapley values and loco. *arXiv:2303.05981*, 2023.

Isabella Verdinelli and Larry Wasserman. Decorrelated variable importance. *Journal of Machine Learning Research*, 25(7):1–27, 2024.

Emmanuel Candes, Yingying Fan, Lucas Janson, and Jinchi Lv. Panning for gold:‘model-x’knockoffs for high dimensional controlled variable selection. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 80(3):551–577, 2018.

Wesley Tansey, Victor Veitch, Haoran Zhang, Raul Rabadan, and David M Blei. The holdout randomization test for feature selection in black box models. *Journal of Computational and Graphical Statistics*, 31(1):151–162, 2022.
Thank you for the comprehensive response! A few follow-up questions:
- For W2, could the authors provide more discussion one when to consider this class of methods over the others added to the related works discussion? 
- For your response Q1, you mentioned, "Under $H_1$, it (the variance estimator) underestimates the actual unknown variance, which is advantageous for power enhancement as demonstrated." Could you explain why that is and where it is demonstrated?