Task-Agnostic Machine-Learning-Assisted Inference

Authors: jiacheng  miao
qiongshi  lu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 18:44:59

Abstract: Machine learning (ML) is playing an increasingly important role in scientific research. In conjunction with classical statistical approaches, ML-assisted analytical strategies have shown great promise in accelerating research findings. This has also opened a whole field of methodological research focusing on integrative approaches that leverage both ML and statistics to tackle data science challenges. One type of study that has quickly gained popularity employs ML to predict unobserved outcomes in massive samples, and then uses predicted outcomes in downstream statistical inference. However, existing methods designed to ensure the validity of this type of post-prediction inference are limited to very basic tasks such as linear regression analysis. This is because any extension of these approaches to new, more sophisticated statistical tasks requires task-specific algebraic derivations and software implementations, which ignores the massive library of existing software tools already developed for the same scientific problem given observed data. This severely constrains the scope of application for post-prediction inference. To address this challenge, we introduce a novel statistical framework named PSPS for task-agnostic ML-assisted inference. It provides a post-prediction inference solution that can be easily plugged into almost any established data analysis routines. It delivers valid and efficient inference that is robust to arbitrary choice of ML model, allowing nearly all existing statistical frameworks to be incorporated into the analysis of ML-predicted data. Through extensive experiments, we showcase our method’s validity, versatility, and superiority compared to existing approaches. Our software is available at https://github.com/qlu-lab/psps.

Main paper: 
 introduction: 
 Leveraging machine learning (ML) techniques to enhance and accelerate research has become increasingly popular in many scientific disciplines [44] . For example, sophisticated deep learning models have achieved remarkable success in predicting protein structure and interactions, which has the potential to significantly speed up the research process, save costs, and revolutionize the field of structural biology [1, 2, 25] . However, recent studies have pointed out that statistical inference using ML-predicted outcomes may lead to invalid scientific discoveries due to the lack of consideration of ML prediction uncertainty in traditional statistical approaches. To address this, researchers have introduced methods that couple extensive ML predictions with limited gold-standard data to ensure the validity of ML-assisted statistical inference [3, 35, 46] . Despite these advances, current ML-assisted inference methods can only address very basic statistical tasks, including mean estimation, quantile estimation, and linear and logistic regression [3, 35] . While the same mathematical principle behind existing ML-assisted inference methods can be generalized to a broader class of M-estimation problems, specific algebraic derivations and computational implementations are required for each new statistical task. Moreover, many tasks, such as the Wilcoxon rank-sum test, do not fit into the M-estimation framework. These issues pose significant challenges to the broad application of ML-assisted inference across various scientific domains.• Classical statistical methods apply algorithm A to only labeled data L = (X L , Y L ), and returns the estimator and its estimated variance [ θ L , Var( θ L )]. Valid confidence intervals and hypothesis tests can then be constructed using the asymptotic distribution of the estimator. However, it ignores the unlabeled data and ML prediction.• Imputation-based methods treat ML prediction f in the unlabeled data as the observed outcome, and apply algorithm A to U = (X U , f U ). We denote the estimator and estimated variance as [ η U , Var( η U )]. This has been shown to give invalid inference results and false scientific findings [3, 35, 36, 46 ]. • ML-assisted inference methods use both L and U as input. These approaches add a debiasing term in the loss function (or estimating equation) for M-estimation problems, thus removing the bias from the imputation-based estimators and producing results that are statistically valid and universally more powerful compared to classical methods [4, 35, 36] .Next, we use an example to provide intuition on ML-assisted inference and our protocol. 
building the intuition with mean estimation: 
 We consider the mean estimation problem, whereθ * = E[Y i ] ≡ arg min θ E[ 1 2 (Y i -θ) 2 ]. The classical method only takes the labeled data Y L as input and yields an unbiased and consistent estimator for θ * : θ L = arg min θ1 n n i=1 1 2 (Y i -θ) 2 = 1 n n i=1Y i . The imputation-based method only takes the unlabeled data f U as input and returns η U = arg min θ1 N n+N i=n+1 1 2 ( f i -θ) 2 = 1 N n+N i=n+1 f i .It is a biased and inconsistent estimator for E[Y i ] if the ML model f is mis-specified.To address this, ML-assisted estimator takes both labeled data (Y L , f L ) and unlabeled data f U as input and adds a debiasing term to the loss function to rectify the bias caused by ML imputation [3, 20, 35, 36] :θMLA = arg min θ 1 2 { ω0 1 N n+N i=n+1 ( fi -θ) 2 -[ ω0 1 n n i=1 ( fi -θ) 2 - 1 n n i=1 (Yi -θ) 2 ] Debiasing term } = ω0 1 N n+N i=n+1 fi -[ ω0 1 n n i=1 fi - 1 n n i=1 yi]Debiasing term, where the modified loss ensures the consistency of the ML-assisted estimator and the weightω 0 = Cov l [Y, f ]/n Var l [ f ]/n+ Varu[ f ]/Nensures that ML-assisted estimator is no less efficient than the classical estimator with arbitrary ML predictions:Var( θ MLA ) = Var( θ L ) - Cov[Y, f ] n Var[ f ]+n 2 Var[ f ]/N ≤ Var( θ L ).Our proposed method is motivated by the observation that the sufficient statistics of the ML-assisted estimator θ MLA and its estimated variance Var( θ MLA ) are the following summary statistics:θ ss = ( 1 n n i=1 y i , 1 n n i=1 f i , 1 N n+N i=n+1 f i ) and Var( θ ss ) =   Var l [Y ]/n Cov l [Y, f ]/n 0 Cov l [Y, f ]/n Var l [ f ]/n 0 0 0 Var u [ f ]/N  Moreover, they can be easily obtained by applying the same algorithm A (mean estimation) to• labeled data with observed outcomeA(Y L ) → [ θ L , Var( θ L )] = [ 1 n n i=1 y i , Var l [Y ]/n]• labeled data with predicted outcomeA( f L ) → [ η L , Var( η L )] = [ 1 n n i=1 f i , Var l [ f ]/n]• unlabeled data with predicted outcomeA(f U ) → [ η U , Var( η U )] = [ 1 N n+N i=n+1 f i , Var u [ f ]/N ]• bootstrap of labeled data A[(Y L , f L ) q , q = 1, . . . , Q] for estimation of Cov( θ L , η L ) = Cov l [Y, f ]/n. Here, (Y L , f L ) q represents the q-th bootstrap of labeled data.Combining these summary statistics for one-step debiasing ω 0 η U -( ω 0 η L -θ L ) recovers θ MLA .To summarize, an algorithm for mean estimation, coupled with resampling, is sufficient for MLassisted mean estimation. This observation inspired us to generalize this protocol for a broad range of tasks. Our protocol illustrated in Fig. 1 only requires three steps: 1) using a pre-trained ML model to predict outcomes for labeled and unlabeled data, 2) applying existing analysis routines to generate summary statistics, and 3) using these statistics in a debiasing procedure to produce statistically valid results in ML-assisted inference.Labeled dataOne-step debiasingStep1: Use pre-trained ML model to predict the gold-standard outcome in both labeled and unlabeled dataApply the analysis routines to obtain the summary statisticsEmploy one-step debiasing to the summary statistics 𝑿 𝓛Figure 1 : Workflow of PSPS for Task-Agnostic ML-Assisted Inference.Our work is closely related to recent methods developed in the literature of ML-assisted inference [3, 4, 20, 35, 37, 46, 56] , and is also related to methods for handling missing data [40, 42] and semi-supervised inference [6, 16, 50, 52] . While current ML-assisted inference methods modify the loss function or the estimating equation, our protocol works directly on the summary statistics. For simple problems such as mean estimation, current methods yield a closed-form solution to the optimization problem. However, for more general statistical tasks, there is no such closed-form solution. Current methods typically require the algebraic form of the loss function, its first-and second-order derivatives, and the variance for the estimator, as well as a newly implemented optimization algorithm to obtain the estimator. We use the logistic regression problem as an example. Here, θ * = arg min θ E[-Y (θX) T -ψ(Xθ)] andψ(t) = 1/(1 + exp(-t)). The ML-assisted estimator is θ MLA = arg min θ 1 N n+N i=n+1 ω[-f i θ T X T i - ψ(X i θ)]-{ 1 n n i=1 ω[-f i θ T X T i -ψ(X i θ)]-1 n n i=1 [-Y i θ T X T i -ψ(X i θ)]} with estimated asymptotic variance A -1 V(ω) A -1 , where A = 1 N +n ( n i=1 ψ ′′ (X i θ) X T i X i + n+N i=n+1 ψ ′′ (X i θ) X T i X i ), V(ω) = n N [ ω 2 Var n (ψ ′ (X i θ) -f i )X T i + Cov N +n (1 -ω)ψ ′ (X i θ)X T i + ( ω f i -Y i )X T i ,and ω needs to be obtained by optimization to minimize the asymptotic variance. In contrast, our protocol simply applies logistic regression A to• labeled data with observed outcomes (X L , Y L ) to obtain [ θ L , Var( θ L )] • labeled data with predicted outcomes (X L , f L ) to obtain [ η L , Var( η L )] • unlabeled data with predicted outcomes (X U , f U ) to obtain [ η U , Var( η U )] • bootstrap of labeled data (X L , Y L , f L ) q , q = 1, . . . , Q for Cov( θ L , η L ),and returns ω T 0 η U -(ω T 0 η L -θ L ), where ω 0 = ( Var( η L )+ Var( η U )) -1 Cov( θ L , η L ).For each new statistical task, as long as an existing analysis routine can produce an estimator that is asymptotically normally distributed, our protocol can be similarly applied. Additionally, the current mathematical principles guiding ML-assisted inference apply solely to M-estimation [3, 4, 20, 35, 56] . Our protocol extends beyond this limitation, addressing all estimation problems with an asymptotically normally distributed estimator.Inference relying solely on summary statistics is widely used in the statistical genetics literature for practical reasons. Summary statistics-based methods have been developed for tasks such as variance component inference and genetic risk prediction [11, 12, 34, 39, 53] . In contrast to our work, these applications do not leverage ML predictions, but instead focus on inference using summary statistics obtained from observed outcomes. An exception is a previous study for valid genome-wide association studies (GWAS) on ML-predicted outcome [36] . However, it focused only on linear regression modeling with application to GWAS. The PSPS framework introduced in this paper aims to extend ML-assisted inference to general statistical tasks.Our work is also related to semi-supervised learning, resampling-based inference, zero augmentation, and false discovery rate (FDR) control methods. Our protocol is designed for estimation and statistical inference using both labeled and unlabeled data, addressing a different problem from semi-supervised learning [55] , which primarily focuses on prediction. Our protocol is inspired by the core principle of resampling-based inference, which replaces algebraic derivation with computation [19] . The main difference is that we focus on how to use ML to support inference, whereas resamplingbased inference focuses on bias and variance estimation, and type-I error control. The idea of zero augmentation has been used in augmented inverse propensity weighting estimators [38] and in handling unmeasured confounders [48] and missing data for U-statistics [14] . These estimators do not incorporate ML, which is fundamental to our work. We also adapt techniques from the FDR literature [7] [8] [9] [10] . Our unique contribution is to use ML to support FDR control, thereby increasing its statistical power, in contrast to classical methods that rely solely on labeled data.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: We have discussed the limitations and future direction of our paper in Conclusion section of the paper.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create 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. 
general protocol for psps: 
 Building on Section 2, we formalized our protocol in Fig. 1 for ML-assisted inference:Algorithm 1 PSPS for ML-assisted inferenceInput: A pre-trained ML model f , labeled data L = (X L , Y L ), unlabeled data U = X U 1:Use the ML model f to predict the outcome in both labeled and unlabeled data.2: Apply the algorithm A in the analysis routine to • labeled data (X L , Y L ) and obtain [ θ L , Var( θ L )]• labeled data (X L , f L ) and obtain [ η L , Var( η L )]• unlabeled data with (X U , f U ) and obtain [ η U , Var( η U )]• Q bootstrap of labeled data (X L , Y L , f L ) q , q = 1, . . . , Q and obtain Cov( θ L , η L ). 3: Employ one-step debiasing to the summary statistics in step2:θ PSPS = ω T 0 η U -( ω T 0 η L -θ L ), where ω 0 = [ Var( η L ) + Var( η U )] -1 Cov( θ L , η L ) and Var( θ PSPS ) = Var( θ L ) - Cov( θ L , η L ) T [ Var( η L ) + Var( η U )] -1 Cov( θ L , η L ) Output: ML-assisted point estimator θ PSPS , standard error Var( θ PSPS ), α-level confidence interval for the k-th coordinate C PSPS α,k = ( θ PSPS k ± z 1-α/2 Var( θ PSPS ) kk ), and (two-sided) p-value 2(1 -Φ(| θPSPS k √ Var( θPSPS) kk ) |)), where Φ is the CDF of the standard normal distribution.The only requirements for our protocol are: i) algorithm A, when applied to labeled data (X L , Y L ), returns a consistent and asymptotically normally distributed estimator of θ * ; ii) labeled and unlabeled data are independent and identically distributed. Under these assumptions, the summary statistics have the following asymptotic properties:EQUATIONwhere η ≡ η(P f ) ∈ R K is defined on (X, f ) ∼ P f , V(•) denotes the asymptotic variance and covariance of a estimator, and ρ = n N . The asymptotic approximation givesV( θ L ) ≈ n Var( θ L ), V( θ L , η L ) ≈ n Cov( θ L , η L ), V( η L ) ≈ n Var( η L ) and V( η U ) ≈ N Var( η U ).Here, we do not require η L and η U to be consistent for θ * , thus allows arbitrary ML model.With the summary statistics following a multivariate normal distribution asymptotically, the debiased estimatorθ PSPS = ω T 0 η U -( ω T 0 η L -θ L )is consistent for θ * and asymptotically normally distributed (Theorem 1). Therefore, by plugging in a consistent estimator for its asymptotic variance V( θ PSPS ) ≈ n Var( θ PSPS ), valid confidence interval and hypothesis testing can be achieved.Remark 1. PSPS is more "task-agnostic" than existing methods in three aspects:1. For M-estimation tasks, currently, only mean and quantile estimation, as well as linear, logistic, and Poisson regression, have been implemented in software tools and are ready for immediate application. For other M-estimation tasks, task-specific derivation of the MLassisted loss functions and asymptotic variance via the central limit theorem are necessary. After that, researchers still need to develop software packages and optimization algorithms to carry out real applications. In contrast, PSPS only requires already implemented algorithms and software designed for classical inference based on labeled data.2. For problems that do not fall under M-estimation but have asymptotically normally distributed estimators, only PSPS can be applied, and all current methods would fail. The principles behind ML-assisted M-estimation do not extend to these tasks.3. Even for M-estimation tasks that have already been implemented, PSPS offers the additional advantage of relying solely on summary statistics. The "task-specific derivations" refer not only to statistical tasks but also to scientific tasks. Real-world data analysis in any scientific discipline often involves conventions and nuisances that require careful consideration. For example, our work is partly motivated by GWAS [43] . Statistically, GWAS is a linear regression that regresses an outcome on many genetic variants. While the regressionbased statistical foundation is simple, conducting a valid GWAS requires accounting for numerous technical issues, such as sample relatedness (i.e., study participants may be genetically related) and population structure (i.e., unrelated individuals of the same ancestry are both genetically and phenotypically similar, creating confounded associations in GWAS). Sophisticated algorithms and software have been developed to address these complex issues [31] . It will be very challenging if all these important features need to be reimplemented in an ML-assisted GWAS framework. With our PSPS protocol, researchers can utilize existing tools that are highly optimized for genetic applications to perform ML-assisted GWAS. This adaptability is not just limited to GWAS, but is a major feature of our approach across scientific domains. PSPS enables researchers to conduct ML-assisted inference using well-established data analysis routines.Remark 2. The "federated data analysis" feature of PSPS refers to the fact that we only require summary statistics as input for inference, rather than individual-level raw data (features X and label Y ). For example, consider a scenario where labeled data is in one center and unlabeled data is in another, yet researchers cannot access individual-level data from both centers simultaneously. Under such conditions, current ML-assisted inference, which relies on accessing both labeled and unlabeled data to minimize a joint loss function, is not feasible. However, PSPS circumvents this issue by aggregating summary statistics from multiple centers, thereby performing statistical inference while upholding the privacy of individual-level data.In this section, we examine the theoretical properties of PSPS. In what follows, P → denotes convergence in probability and D → denotes convergence in distribution. All proofs are deferred to the Appendix A.The first result shows that our proposed estimator is consistent, asymptotically normally distributed, and uniformly better in terms of element-wise asymptotic variance compared with the classical estimator based on labeled data only.Theorem 1. Assuming equation (1) holds, then θ PSPS P → θ * , andn 1/2 ( θ PSPS -θ * ) D → N 0, V( θ PSPS ) , where V( θ PSPS ) = V( θ L )-V( θ L , η L ) T (V( η L )+ρV( η U )) -1 V( θ L , η L ). Assume the k-th column of V( θ L , η L ) is not a zero vector and at least one of V( η L ) and V( η U ) are positive definite, then V( θ PSPS ) kk ≤ V( θ L ) kk . With V( θ PSPS ) P → V( θ PSPS ), lim n P(θ * k ∈ C PSPS α,k ) = 1 -α.V( θ PSPS ) can be obtained by applying the algebraic form of V( θ PSPS ) using the bootstrap estimators for V( θ L ), V( η L ), V( θ L , η L ), and V( η U ). The regularity conditions for consistent bootstrap variance estimation are outlined in Theorem 3.10 (i) of [41] . We also refer readers to [21] , which showed that bootstrap-based variance provides valid but potentially conservative inference.This result indicates that a greater reduction in variance for the ML-assisted estimator is associated with larger values of V( θ L , η L ) and smaller values of V( η L ), V( η U ), and ρ. The variance reduction term[V( θ L , η L ) T (V( η L ) + ρV( η U )) -1 V( θ L , η L )] kk can also serve as a metric for selecting the optimal ML model in ML-assisted inference.Our next result shows that three existing methods, i.e., PPI, PPI++, and PSPA, are asymptotically equivalent to PSPS with different weighting matrices. A broader class for consistent estimator of θ* is θ(ω) = ω T η U -(ω T θ L -η L ), where ω is a K × K matrix. The consistency of θ(ω) for θ * only requires ω T ( η U -η L ) P → 0. Since ( η U -η L )P → 0, assigning arbitrarily fixed weights for will satisfy the condition. However, the choice of weights influences the efficiency of the estimator as illustrated in Proposition 2 later. Proposition 1. Assuming equation (1) and regularity condition for the asymptotic normality of current ML-assisted estimator holds. For any M-estimation problem, we haven 1 2 ( θ (diag(ω ele )C)-θ PSPA ) D → 0, n 1 2 ( θ (diag(ω tr )C)-θ PPI++ ) D → 0, n 1 2 ( θ (diag(1)C)-θ PPI ) D → 0.Here, ω ele = [ω ele,1 , . . . , ω ele,K ] T ∈ R K and ω ele,k minimizing the k-th diagonal element of V( θ(ω)), ω tr is a scalar used to minimize the trace of V( θ(ω)), and C is a matrix associated with the second derivatives of the loss function in M-estimation, with further details deferred to Appendix A.This demonstrates that for M-estimation problems, our method is asymptotically equivalent to PSPA, PPI++, and PPI with the respective weights diag(ω ele )C, diag(ω tr )C, and diag(1)C. Therefore, PSPS can be viewed as a generalization of these existing methods.Our third result shows that the weights used in the Proposition 1 are not optimal. Instead, our choice of ω 0 represents the optimal smooth combination of ( θ L , η L , η U ) in terms of minimizing the asymptotic variance, while still preserving consistency.Proposition 2. Suppose n 1/2 (g( θ L , η L , η U ) -θ * ) D → N (0, Σ g ) and g is a smooth function, then Σ g kk ≥ Σ PSPS kkTogether with Proposition 1, our results demonstrate that our protocol provides a more efficient estimator compared to existing methods for the M-estimation problems. Furthermore, the applicability of our protocol is not limited to M-estimation and only requires summary statistics as input. It also indicates that in a setting of federated data analysis [24] where individual-level data are not available, PSPS proves to be the optimal approach for combining shared summary statistics. Remark 3. PPI++ [4] employs a power-tuning scalar for variance reduction in ML-assisted inference. This scalar is obtained by minimizing the trace or possibly other scalarization of the estimator's variance-covariance matrix. However, the asymptotic variance of PSPS is always equal to or smaller than that of PPI++, irrespective of the scalarization chosen by researchers. This advantage arises because PSPS utilizes a K × K power tuning matrix, ω, for variance reduction, where K represents the dimensionality of parameters. This matrix facilitates information sharing across different parameter coordinates, thereby enhancing estimation precision. The choice of weighting matrix in PSPS also allows for element-wise variance reduction, reducing each diagonal element of the variance-covariance matrix. In contrast, the single scalar in PPI++ can only target overall trace reduction or variance reduction of a specific element. A detailed example is provided in Appendix B. Only in one-dimensional parameter estimation tasks, such as mean estimation, PPI++ and PSPS exhibit the same asymptotic variance.We also provide several extensions to ensure the broad applicability of our method.Here, we relax the assumption that the labeled data and unlabeled data are independent. When they are not independent, this can lead to the non-zero covariance between the η L and η U . Consider a broader class of summary statistics asymptotically satisfyingn 1/2   θ L -θ * η L -η η U -η   D → N    0 K 0 K 0 K ,     V( θ L ) V( θ L , η L ) V( η L , η U ) V( θ L , η L ) V( η L ) √ ρV( θ L , η U ) V( η L , η U ) √ ρV( θ L , η U ) ρV( η U )       We can similarly employ the one-step debiasingθ no-indep PSPS = ω T 0 η U -ω 0 ( η L -θ) where ω 0 = ( V( θ L , η L ) -V( η L , η U )) T ( V( η L ) + V( η U ) -2 V( θ L , η U )) -1 and Var( θ no-indep PSPS ) = Var( θ L ) - ( V( θ L , η L ) -V( η L , η U )) T (V( η L ) + V( η U ) -2 V( θ L , η U )) -1 ( V( θ L , η L ) -V( η L , η U )).The theoretical guarantees of the proposed estimator can be similarly derived by Theorem 1.The other assumption of our approach is that the labeled and unlabeled data are identically distributed so that we can ensure η L -η U P → 0 and validity of PSPS results. To address the potential violation of this assumption, we introduce a sensitivity analysis with hypothesis testing for the null H 0 :η L,k = η U ,k with test statistics η L,k -η U ,k √ Var( η L,k )+ Var( η U ,k ) D → N (0, 1) to assess if η L,k and η U ,k are significantly different.Here, the subscript k indicates the k-th coordinate. We recommend using p-value < 0.1 as evidence for heterogeneity and caution the interpretation of results from ML-assisted inference.The output from PSPS can be used for ML-assisted FDR control, achieving greater statistical power compared to classical FDR control methods that solely rely on labeled data. We refer to our approach as PSPS-BH and PSPS-knockoff. Briefly, PSPS-BH processes the p-value from ML-assisted linear regression through the Benjamini-Hochberg (BH) procedure [9] , while PSPS-knockoff utilizes the ML-assisted debiased Lasso coefficient [23, 51] in the ranking algorithm of knockoff [7] . We present our algorithm in Appendix C and evaluate their performance using experiments in Section 4.We have discussed ML-assisted inference with outcomes predicted by ML models. Here, we note that PSPS can also be applied when either features alone are predicted or both features and outcomes are predicted. The key idea is that the difference between point estimators obtained from applying A to predicted features in both labeled and unlabeled datasets is a consistent estimator for zero. This enables zero augmentation for estimators from observed features and outcomes. To implement this, modify step 2 in Algorithm 1 to apply A to predicted features in both labeled and unlabeled data. A similar approach is applicable when both features and outcomes are predicted.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?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 have provided the code in the Supplementary Material.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
numerical experiments and real data application 4.1 simulations: 
 We conduct simulations to assess the finite sample performance of our method. Our objectives are to demonstrate that 1) PSPS achieves narrower confidence intervals when applied to statistical tasks already implemented in existing ML-assisted methods; 2) when applied to statistical tasks that have not been implemented for ML-assisted inference, PSPS provides confidence intervals with narrower width but correct coverage (indicating higher statistical power) compared to classical approaches rely solely on labeled data; 3) PSPS provides well-calibrated FDR control and achieves higher power compared to classical methods only using labeled data.Tasks that have been implemented for ML-assisted inference We compare PSPS with the classical method using only labeled data, the imputation-based method using only unlabeled data, and three ML-assisted inference methods PPI, PPI++, and PSPA [3, 4, 35] on mean estimation and linear and logistic regression. We defer the detailed data-generating process to Appendix D. In short, we generated outcome Y i from feature X 1i and X 2i , and obtained the pre-trained random forest that predict Y i using X 1i and X 2i . We have 500 labeled samples (X 1i , Y i , f i ) , and unlabeled samples (X 1i , f i ) ranged from 1,000 to 10,000. Our goal is to estimate the mean of Y i , as well as the linear and logistic regression coefficient between Y i and X 1i . Fig. 2a-c show confidence interval coverage and Fig. 2d-f show confidence interval width. We find that the imputation-based method fails to obtain correct coverage, while all others including PSPS have the correct coverage. PSPS has narrower confidence intervals compared to the classical method and other approaches for ML-assisted inference. Tasks that have not been implemented for ML-assisted inference Next, we consider several commonly used statistical tasks that currently lack implementations for ML-assisted inference including quantile regression [27] , instrumental variable (IV) regression [5] , negative binomial (NB) regression [22] , debiased Lasso [51] , and the Wilcoxon rank-sum test [28] . Similar to our previous simulations, we utilize labeled data, unlabeled data, and a pre-trained ML model. Detailed simulation settings are deferred to the Appendix D. Our goal is to estimate the regression coefficient between Y i and X 1i for quantile (at quantile level 0.75), IV, and NB regression, between Y i and high dimensional features X i ∈ R 150 for debiased Lasso, and to perform hypothesis testing on the medians of two independent samples Y i |X 1i = 1 and Y i |X 1i = 0 using the Wilcoxon rank-sum test.Fig. 3a-d show confidence interval coverage and Fig. 3f -i show confidence interval width for parameter estimation. Fig. 3e and Fig. 3j show the type-I error and statistical power for the Wilcoxon rank-sum test. We found that the imputation-based method fails to obtain correct confidence interval coverage and shows inflated type-I error, while PSPS and classical method have the correct coverage and wellcalibrated type-I error control. PSPS has narrower confidence intervals width in all settings, and higher statistical power for the Wilcoxon rank-sum test compared to classical methods. Confidence intervals become narrower as unlabeled sample size increases, indicating higher efficiency gain.We evaluate the finite sample performance of PSPS-BH and PSPS-knockoff in controlling the FDR compared with classical and imputation-based methods as baselines. We consider low-dimensional(K < n) and high-dimensional(K > n) linear regressions for PSPS-BH and PSPS-knockoff, respectively. We simulate the data such that only a proportion of the features are truly associated with the outcome. The data generating process is deferred to Appendix D. Our goal is to select the associated features while maintaining the target FDR level. Imputation-based method failed to control FDR in either low-dimensional or high-dimensional settings. Classical approach, PSPS-BH, and PSPS-knockoff effectively controlled in both lowdimensional and high-dimensional settings. PSPS-BH, and PSPS-knockoff achieve higher statistical power compared to the classical method.These simulations demonstrate that PSPS outperforms existing methods and can be easily adapted for various statistical tasks not yet implemented in current ML-assisted inference methods.We employed our method to carry out ML-assisted quantile regression to identify genetic variants associated with the outcome variability (vQTL) of bone mineral density derived from dual-energy X-ray absorptiometry imaging (DXA-BMD) [33] . DXA-BMD is the primary diagnostic marker for osteoporosis and fracture risk [15, 54] . Identifying vQTL for DXA-BMD can provide insights into the biological mechanisms underlying outcome plasticity and reveal candidate genetic variants involved in potential gene-gene and gene-environment interactions [29, 32, 45, 47, 49] . We focused on total body DXA-BMD, which integrates measurements from multiple skeletal sites. We used data from the UK Biobank [13] , which includes 36,971 labeled and 319,548 unlabeled samples with 9,450,880 genetic variants after quality control. We predicted DXA-BMD values in both labeled and unlabeled samples using SoftImpute [30] with 466 other variables measured in the UK Biobank. Prediction in the labeled sample was implemented through cross-validation to avoid overfitting. The implementation detail is deferred to Appendix D. We used the BH procedure to correct for multiple testing and considered FDR < 0.05 as the significance threshold.No genetic variants reached statistical significance under the classical method with only labeled data. PSPS identified 108 significant variants with FDR < 0.05 spanning 5 independent loci, showcasing the superior statistical power of PSPS (Fig. E.2 and Table E .1). Notably, these significant vQTL cannot be identified by linear regression [36] , indicating the different genetic mechanisms controlling outcome levels and variability for DXA-BMD.We compared the computational efficiency of PSPS with existing methods using a dataset of 500 labeled and 10,000 unlabeled data points. Results are shown in Table E .2. While PSPS is slower due to resampling, its overall runtime is still relatively short. 
conclusion: 
 We introduced a simple, task-agnostic protocol for ML-assisted inference, with applications across a broad range of statistical tasks. We established the consistency and asymptotic normality of the proposed estimator. We further introduced several extensions to expand the scope of our approach.Through extensive experiments, we demonstrated the superior performance and broad applicability of our method across diverse tasks. Our protocol involves initially generating summary statistics using computationally efficient software tools in scientific data analysis, followed by integration of summary statistics to produce ML-assisted inference results, which achieves high computational efficiency while maintaining statistical validity. Future work could focus on developing a fast resampling algorithm to further improve computational efficiency. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: The detail of the experimental can be found in Appendix D.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We reported the confidence interval width and p-value for statistical significance throughout the paper.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]Justification: The detail of the experiments compute resources can be found in Appendix D.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: We have followed the code Of ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes]Justification: We have discussed the positive societal impacts of our method is accelerating scientific in the abstract and introduction.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: Our paper introduced a statistical method and therefore 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: We have cited all paper that produced the code package or 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. 
 
 Appendix: A.1 Proof the Theorem 1Proof. SinceEQUATIONθ L P → θ * and η U -η L P → 0. Given weights ω 0 = ( V( η L ) + ρ V( η U )) -1 V( θ L , η L ), where the variances are consistently estimated, Slutsky's theorem implies that, ω T 0 ( η U -η L )P → 0.Also by Slutsky's theorem,EQUATIONwhich completes the proof of consistency.By multivariate delta methods, denoting h([x, y, z] T ) = x + ω T 0 (z -y), we have ∇h([x, y, z] T ) = [1, -ω 0 , ω 0 ], therefore by the consistency of ω 0 ,EQUATIONEQUATIONwhich completes the proof of asymptotic normality.EQUATIONHere, by our assumption, (V( η L ) + ρV( η U )) -1 is a positive definite matrix. Therefore, quadratic form V( θ L , η L ) T :,k (V( η L ) + ρV( η U )) -1 V( θ L , η L ) :,k is non-negative and is zero if only all elements of V( θ L , η L ) :,k is zero, which completes the proof of element-wise variance reduction.Proof. GivenEQUATIONthe asymptotic variance of θ(ωEQUATIONWe first define the M-estimation (Z-estimation) problem. The goal is to estimate a K-dimensional parameter θ * defined through an estimating equationEQUATIONwhere ψ(•, •; θ) is a user-defined function. By the theory of Z-estimation and a recent paper on MLassisted inference [35] , we haveV( θ L ) = A -1 θ * M 1 A -1 θ * , V( θ L , η L ) = A -1 η M 4 A -1 θ * , V( η L ) = A -1 η M 2 A -1 η , and V( η U ) = A -1 η M 3 A -1 η .Here,M 1 = Var l [ψ(Y, X; θ * )], M 2 = Var l [ψ( f , X; θ * )], M 3 = Var u [ψ( f , X; θ * )], M 4 = Cov l [ψ(Y, X; θ * ), ψ( f , X; θ * )], A θ * = E[∂ψ(Y, X; θ * )/∂θ], A η = E[∂ψ( f , X; η)/∂η], and ρ = n N .Rewritten V( θ(ω)) using the above notation, we haveEQUATIONPlug inEQUATIONinto the equation above, we haveEQUATION)To connect our protocol with existing methods, we define EQUATIONEQUATIONBy the theory of PSPA, PPI++, and PPI paper [3, 4, 35] , we haveEQUATION)EQUATION)EQUATIONBased on the proof of Theorem 1, we have theEQUATION)EQUATION)EQUATIONDenote C = A -1 θ * A η , we have PSPA, PPI++, and PPI is asymptotically equivalent to diag(ω * ele )C, diag(ω * tr )C, and diag(1)C, respectively. This completes the proof.Proof. We apply the first-order Taylor expansion to g( θ L , η L , η U ) around (θ * , η, η):EQUATIONwhere we used ∇ η,1 and ∇ η,2 to denote the gradient of g(θ * , η, η) w.r.t the first and second η, respectively.This can be written as a linear function of (θ * , η, η):EQUATIONwhereµ = g(θ * , η, η) -∇ θ * g(θ * , η, η)θ * -2∇ η g(θ * , η, η)η, β 1 = ∇ θ * g(θ * , η, η), β 2 = ∇ η,1 g(θ * , η, η), β 3 = ∇ η,2 g(θ * , η, η).Since we require g( θ L , η L , η U ) P → θ * , we have µ = 0, β 1 = 1, and β 2 + β 3 = 0. This leads toEQUATIONGivenEQUATIONwe haveEQUATIONwhich is a quadratic function of β 3 and achieves it minimum whenβ 3 = (V( η L ) + ρV( η U )) -1 V( θ L , η L ) = ω 0 .This completes the proof.Consider a linear regression with two predictors: Since θPSPS,1 involves two zero-augmentation terms (i.e., -w 1 η1L +w 1 η1U and -ω 12 η2L + ω 12 η2U ), its asymptotic variance should be less than or equal to that of PPI++ with one augmentation term. Therefore, PSPS borrows information from both coordinates, but PPI++ is restricted to information from only the first coordinate. Although the PPI++ can be used under a different scalarization, it still contains one augmentation term.Y ∼ θ 1 X 1 + θ 2 X 2 .Algorithm 2 PSPS-BH for linear regression Input: Labeled data L = (X L , Y L , f L ), unlabeled data U = (X U , f U ), FDR level q ∈ (0,:= p (k) , where k = max k = 1, . . . , K : p (k) ≤ k K q Output: Discoveries S = k : p k ≤ τ BH q Algorithm 3 PSPS-knockoff with debiased Lasso Input: Labeled data L = (X L , Y L , f L ), unlabeled data U = (X U , f U ), FDR level q ∈ (0, 1).1: Obtain the augmented labeled and unlabeled data as L = (X L , XL , Y L , f L ) and Ũ = (X U , XU , f U ) where XL ← knockoff-sample(X L ) and XU ← knockoff-sample(X U ). 2: Calculate the PSPS debiased Lasso coefficient β PSPS-DLasso ← PSPS-DLasso( L, Ũ)3: W k for k = 1, • • • , K ← knockoff-statistic( β PSPS-DLasso ) 4: Set the cutoff τ knockoff q = t > 0 : #{k:W k ≤-t} #{k:W k ≥t}∨1 ≤ q Output: Discoveries S = k : M k > τ knockoff qHere, we employ second-order multivariate Gaussian knockoff variables for knockoff-sample and use the difference between the absolute values of the k-th feature and its knockoff coefficient as the knockoff-statistic. Alternative choices for these two steps can also be readily integrated into our algorithm [26] .Here, we provide the details for our simulation. All our simulation is run in R with version 4.2.1 (2022-06-23) in a MacBook Air with an M1 chip. For all the simulations, the ground truth coefficients are obtained using 5 × 10 4 samples; A pre-trained random forest with 500 trees to grow is obtained from hold-out data. We bootstrap the labeled data for 200 times for covariance estimation. All simulations are repeated 1000 times.• Mean estimation, Linear regression, and Quantile regression: We simulate the data fromY i = X 1i β 1 + X 2i β 2 + X 2 1i β 2 + X 2 2i β 2 + ϵ i, where X i1 and X 2i are independent simulated from N (0, 1),β 1 = √ 0.08, β 2 is set to be the value such that X 2i β 2 + X 2 1i β 2 + X 2 2iβ 2 explains 81% of the outcome variance, and ϵ i is simulated from a mean zero normal distribution with variance such that Var(Y i ) = 1. We use X 1i and X 2i as features to predict the Y i in the random forest. We consider labeled data with 500 individuals, and unlabeled data with sample size 1000, 2500, 5000, or 10000.• Logistic regression: We simulate the data fromY i = 1( Ỹi > median ( Ỹi )), where Ỹi = X 1i β 1 + X 2i β 2 + X 2 1i β 2 + X 2 2i β 2 + ϵ i, where X i1 and X 2i are independent simulated from N (0, 1), β 1 = √ 0.08, β 2 is set to be the value such that X 2i β 2 + X 2 1i β 2 + X 2 2i β 2 explains 81% of the outcome variance, and ϵ i is simulated from a mean zero normal distribution with variance such that Var( Ỹi ) = 1. We use X 1i and X 2i as features to predict the Y i in the random forest. We consider labeled data with 500 individuals, and unlabeled data with sample size 1000, 2500, 5000, or 10000.• Instrumental variable (IV) regression: We simulate the data byEQUATIONEQUATIONEQUATIONEQUATIONWe use X 1i and X 2i as features to predict the Y i in the random forest. We consider labeled data with 500 individuals, and unlabeled data with sample size 1000, 2500, 5000, or 10000.The Z i is used as a instrument for X 1i . • Negative binomial (NB) regression: We simulate the data byEQUATIONEQUATIONY i = NegativeBinomial (s = k, µ = µ i ), where s is the dispersion parameter and µ is the rate.We use X 1i and X 2i as features to predict the Y i in the random forest. We consider labeled data with 500 individuals, and unlabeled data with sample size 1000, 2500, 5000, or 10000.• Debiased Lasso: We simulate the data byEQUATIONθ 1 , . . . , θ 15 = 0.9 √ 15 ; θ 16 , . . . , θ 200 = 0 (42)EQUATIONWe use X 1i , . . . , X 200i as features to predict the Y i in the random forest. We consider labeled data with 100 individuals, and unlabeled data with sample size 1500, 2000, 2500, or 3000.• Wilcoxon rank-sum test: We simulate the data byEQUATIONEQUATIONwhere β 1 = √ 0.01 to power simulation and β 1 = 0 for type-I error simulation, β 2 is set to be the value such that β 2 X 2i + β 2 X 2 2i explains 81% of the outcome variance, and τ ϵ is set to be the value such that Var(Y i ) = 0.We use X 1i and X 2i as features to predict the Y i in the random forest. We consider labeled data with 500 individuals, and unlabeled data with sample size 1000, 2500, 5000, or 10000.• Benjamini-Hochberg (BH) procedure: We set K = 150 generate the features independently from N (0, Σ), where Σ is a symmetric Toeplitz matrix that has the structure:EQUATIONThe correlation r is set to be 0.25. We then simulate the outcome Y i = 150 k=1 X ki β k + ϵ i , where we randomly sample 15 β k to be 0.15 and let all other remaining β k to be 0. ϵ i is simulated from a mean-zero normal distribution with variance set to the value such that Var(Y i ) = 1. We further generate Z i = 0.9Y i + 150 k=1 X ki θ k + γ i , where θ k = 0.15 for all k = 1, . . . , 150. We use Z i as features to predict the Y i in the random forest. We consider labeled data with 500 individuals, and unlabeled data with sample size 5000.• knockoff: We used the same setting as described in the BH procedure above to generate the data. The only difference is that we set β k = 0.5 for features associated with the outcome and considered labeled data consisting of 100 individuals, along with unlabeled data comprising a sample size of 1000.Our prediction pipeline comprises two components: prediction for unlabeled data and prediction for labeled data. To predict bone mineral density in unlabeled data, we first selected predictive features by 1) calculating the correlation of bone mineral density with 466 other variables (sample size > 200,000 from UK Biobank) using labeled data and 2) selecting the top 50 variables with the highest correlations as inputs for the SoftImpute algorithm [30] to predict bone mineral density in the unlabeled data. For the labeled data, we employ a similar approach but incorporate 10-fold cross-validation to prevent overfitting. We select the predictive variables and train the SoftImpute model using 90% of the labeled data. We then perform predictions on the remaining 10% in each fold and repeat this process 10 times across all folds. Justification: We provided extensive theoretical and experimental evidence to support the main claims in our paper.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not 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.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 


Decision: The paper proposes a framework for statistical inference that combines labeled, unlabeled through a post-inference solution that can be used with different data analysis routines. The reviewers unanimously recommend its acceptance. They noted that the proposed method is simple and flexible, the evaluations were comprehensive through theoretical and empirical analysis and highlighted its application to important areas. I therefore recommend its acceptance.

Reviews: review: paper_summary: paper_summary: The paper introduces a task-agnostic approach to inference with machine learning predictions. The basic idea follows a similar recipe as prediction-powered inference and related recent papers, but it makes use of resampling instead of the CLT with a plug-in estimate of the asymptotic covariance to avoid relying on analytical problem-specific derivations and expressions.
strength_weakness: strength_weakness: I agree with your point that the recent methods require task-specific derivations. However, it's worth noting that, since for M-estimators we know we get asymptotic normality, we can use the same estimators but instead of deriving the asymptotic variance through the CLT we can use resampling to estimate the asymptotic variance. This is just to say that for M-estimators we can use the old estimators and get inference without problem-specific derivations. (Your other criticisms still apply.)

In L144 it says that resampling-based inference focuses on bias and variance estimation. I don't really agree with this. Plenty of resampling-based inference focuses on type I error control (confidence intervals and p-values).

The mean-estimation result in Section 2.2 is completely borrowed from prior work. Please make this clear. Otherwise it looks like these are new results in this paper.

In Algorithm 1 it is unclear what the actual output is. Please write it out using mathematical symbols.

Regarding Proposition 2, I would say that the claim looks a bit too strong. In PPI++ the authors state that they give the example with minimizing the trace as an example, but that other scalarizations of the covariance are clearly possible. Your Proposition 2 just chooses a different scalarization for the comparison so it is clearly better than the trace example from PPI++. But the PPI++ argument would give the same asymptotic variance under your scalarization.

It would be helpful to give an example of when the p-values in Proposition 3 will be PRDS.

Some stylistic comments:
- The grammatically correct way to spell the title would be "Task-Agnostic Machine-Learning-Assisted Inference."
- in abstract: constraints -> constrains
- Very often throughout the paper you say "standard error" but you write the variance. For example, see L68 or L73. Please be consistent: either say variance and write Var(...) or write out the standard error.
- In L75, you are clearly describing a particular class of ML-assisted inference methods. Not every method that is ML-assisted follows that description. Add references so it's clear what you are referring to.
- L91: motivated the observation -> motivated by the observation
- L99: inputs -> inputting
- There are other typos and minor stylistic issues. Please go through the paper carefully.
questions: questions: - Can you elaborate on the point about the privacy-preserving feature of your method? I didn't understand the main point. A particular use case would be helpful.
- I'm surprised that PSPS and PPI++ are not getting the same interval widths in linear and logistic regression in Figure 2. Is it because you are tuning PPI++ for the trace objective, as discussed above? The two methods (assuming the right tuning) should have exactly the same asymptotics in this problem. This should also be made clear in words.
- I'm curious about the details behind the real-data application. Can you elaborate on how you computed the predictions and how you used cross-validation to avoid overfitting?
limitations: limitations: I don't think a further discussion is necessary.
paper_summary: paper_summary: The paper proposes a novel statistical framework for ML-assisted inference. It describes how labeled data, together with unlabeled data and a pre-trained ML model, can be used for statistical inference. The paper establishes the asymptotic properties and optimality of the framework and then evaluates it empirically on simulated and real datasets.
strength_weakness: strength_weakness: - The limitations could have been discussed more thoroughly.
questions: questions: - What are the limitations of the proposed framework?
limitations: limitations: The authors state that the limitations are discussed in the Conclusion section of the paper, but it seems they are not discussed there.
paper_summary: paper_summary: This paper proposes a new unified framework for ML-assisted inference that reduces the general problem to essentially one of estimating normal means, and then applies simple operations (and a bootstrap step) to solve the normal means problem. In addition to the unifying framework’s simplicity, a key result is that the asymptotic efficiency of the proposed estimator dominates that of existing works under mild conditions.
strength_weakness: strength_weakness: 1. Theorem 1 doesn’t give asymptotically valid inference, since it doesn’t address variance estimation. If the asymptotic variance were known, this result would be sufficient, but it is not, and what is used instead is an estimate of the variance, which must be proved to be consistent in order to conclude asymptotic validity of Algorithm 1’s confidence interval and p-value. This oversight is relevant in other areas, such as eq (35) in the appendix is unjustified and does not follow from eq (34).
2. Prop 1 as stated does not say that these methods all have the same asymptotic variance, which seems to be the implication of the rest of section 3.2 after Prop 1. It seems the proof does prove the right thing, it’s just the statement is wrong (it is missing \sqrt{n} multiplying the LHS of each of the three limits).
3. Proof of prop 4 seems to just be rehashing the proof of knockoffs, except for a single sentence in lines 531-532 which simply states, without proof, the most important part, and the only part that involves the proposed method.

In my opinion, the paper cannot be published without these issues addressed, hence my current score. Item 2 seems to just be a typo, but item 1 is absolutely critical and central, and item 3 is just not rigorous, so it should either be fixed or deleted (personally, I don’t think section 4 strengthens the paper much, since it’s just applying known FDR control techniques to the output of the proposed method, so the validity should just follow from the validity of the proposed method and that of the known FDR control techniques). If they are addressed (and no other reviewers raise other issues that cause me concern), my score would go up quite a bit, as I think this paper has significant strengths.
questions: questions: In what sense is your method more task-agnostic than others? It seems to still require \mathcal{A}, which is task-specific. The authors mention in lines 129-131 that works have been developed that are general to M-estimators, but many (all?) the simulation tasks in Figure 3 which claim to be ones which “have not been implemented for ML-assisted inference” are M-estimation tasks (or reduce to them), so why couldn’t other methods have been applied to them? Explaining this more clearly would raise my score.
limitations: limitations: The authors have not addressed computational efficiency, which is a critical property—how fast is it relative to the other methods under consideration? Addressing this would raise my score.
paper_summary: paper_summary: This paper proposes a post-prediction inference solution (PSPS) that can be adapted into various established data analysis routine and delivers valid and efficient inference for most ML models. In particular, the paper uses both labelled and unlabelled data to derive estimators that are consistent and efficient and only utilizes 1st and 2nd order summary statistics.
strength_weakness: strength_weakness: - The paper claims the proposed method offers better statistical power in various statistical tasks from mean estimation to quantile regression through various numerical experiments in Section 5.1. However, it is unclear how much of the statistical power is offered by the nature of transductive inference and how much is offered by the proposed protocol. For example, [1] shows that incorporating unlabelled data achieves lower test error than purely labelled data already. Could the author demonstrate the unique advantages offers by the proposed method?
- The protocol proposed relies on assumptions on 1) i.i.d. distributions between labelled and unlabelled data and 2) finding an algorithm applied to labeled data returns a consistent and asymptotic normally distributed estimators. I am wondering how easy is it to find an algorithm that can produce such estimators? 
- The authors list one of the key features of the method is privacy-preserving, but did not discuss this point in the remaining paper. 

[1] Chapelle, O., Vapnik, V., & Weston, J. (1999). Transductive inference for estimating values of functions. Advances in Neural Information Processing Systems, 12.
questions: questions: Please see above.
limitations: limitations: Yes.

Comments: comment: Thank you for addressing the concern on the proposed method's applicability - I am raising the score accordingly while maintaining low confidence.
Thank you, this does address my main soundness concerns, and I am raising my score accordingly.
Thanks for your reply!

W2: We would like to clarify that our framework does not require the consistent and asymptotically normal estimators for ML algorithms that are commonly expected in statistics. In fact, our approach can accommodate any "black box" ML algorithm. In our setting, the ML algorithm is used to impute (predict) labels in unlabeled data. These predictions are then used as input to a statistical method (algorithm), such as linear regression, to solve a statistical problem, such as estimating the effect of DNA on height. The requirement for consistent and asymptotically normal estimators applies only to the statistical algorithm, not to the ML algorithm used for label prediction.

We have updated our paper to explicitly state that we did not impose any constraints on the ML algorithms used for label prediction.
Thank you for your valuable comments!

**W1:** We have included in the manuscript the formal regularity conditions required for consistent bootstrap variance estimation. These conditions are detailed in Theorem 3.10 (i) from `Shao, Jun, and Dongsheng Tu. The jackknife and bootstrap. Springer Science & Business Media, 2012.`, which proves the consistency of the bootstrap variance estimator.

Below is the detailed theorem:

We assume that $X _1, \cdots, X _n$ are i.i.d random p-dimensional vectors from distribution $F$. Let $T _n = T _n(X _1, \cdots, X _n)$ be a n estimator of an unknown parameter $\theta$, and $\Re _n = \sqrt{n}(T _n - \theta) \sim N(0, \sigma_n^2)$. Let $\{X _1^*, \ldots, X _n^*\}$ be a bootstrap sample from the empirical distribution $F _n$ based on $X _1, \cdots, X _n$, $T_n^*=T_n\left(X_1^*, \ldots, X_n^*\right)$ and $\Re_n^*=\sqrt{n}\left(T_n^*-T_n\right)$. Denote the bootstrap variance estimator $v _{\text {boot}} = \text{Var}(T_n^*)$.

Theorem 3.10 (i) Let $T _n=T\left(F _n\right)$, assume that $T$ is $\rho _{\infty}$-Fréchet differentiable at $F$, and 
$
\max _{i_1, \ldots, i _n}\left|T _n\left(X _{i _1}, \ldots, X _{i _n}\right)-T _n\right| / \tau _n \rightarrow _{a . s .} 0
$
where the maximum is taken over all integers $i _1, \ldots, i _n$ satisfying $1 \leq$ $i _1 \leq \cdots \leq i _n \leq n$, and $\{ \tau _n \}$ is a sequence of positive numbers satisfying $\liminf _n \tau _n>0$ and $\tau _n=O\left(e^{n^q}\right)$ with a $q \in (0, \frac{1}{2})$., then $v _{\mathrm{boot}} / \sigma _n^2 \rightarrow _p 1$, where $\sigma _n^2=n^{-1} E[\phi _F(X _1)]^2>0$

We have also added a remark in our paper to refer the reader to the theoretical results in a more recent paper `Hahn, Jinyong, and Zhipeng Liao. "Bootstrap standard error estimates and inference." Econometrica 89.4 (2021): 1963-1977.` for bootstrap variance estimation. This paper (THEOREM 1) shows that common bootstrap based standard error in fact leads to a valid but (potentially conservative) inference.

**W3:** Thank you for bringing this to our attention. Since we have decided to remove Section 4 from the paper and keep only the remarks in Section 3, we will no longer include the theoretical results related to this section in the paper. Instead, we will empirically verify the performance of ML-assisted FDR control (as we have done in Section 5 of the paper). As the reviewer pointed out, this does not affect the main contribution of our paper, which is the compelling feature of "task-agnostic" ML-assisted inference.

**Task-agnostic feature of PSPS:**
In Section 3.1, we have added a paragraph specifically highlighting why PSPS is considered "task-agnostic." This addition includes the three bullet points from our previous rebuttal that clearly delineate the relevant scenarios.

**Runtime comparison:**
We have added the runtime comparison into the Section 5 (Numerical experiments and real data application) and the relevant discussions in the Section 6 (Conclusion) in the paper.

We hope our responses address your concerns, and please do not hesitate to let us know if you have any further questions.
W1: Thank you for the clarification. I agree that the current work addresses different problems with [1]. 

W2: I agree consistent and asymptotically normal estimators are common in statistical sciences but not convinced it is common for current ML algorithms, e.g., LLM, given the paper focuses on ML-assisted inference. 

W3: Thank you for the clarification. 

In general, I agree that this paper proposes a new solution that improves statistical consistency and correct confidence coverage using unlabeled data, however I am not certain the impact of contribution given I am not an expert in this area. Therefore I will raise my score but lower my confidence to reflect this.
Thank you for the response! It was very helpful and clarifying.
W1: I agree that IF the variance estimator is consistent, then everything works out easily via Slutsky's--this is not the issue. The main issue is that consistency of the variance estimator is not proved, including in the rebuttal, which simply says "bootstrap estimators of variance is consistent under certain regularity conditions". This is not a proof--what are these conditions, and what papers prove that consistency under those conditions?

W2: Thank you.

W3: Consistency to 0 of the debiased lasso does not imply asymptotic symmetry of the null knockoff statistics. E.g., a Gaussian with mean 1/n and standard deviation 1/n converges to 0 but is not asymptotically symmetric about 0.

As these soundness issues (primarily W1+W3) were the main reason for my low score, I am not revising my score at this time. I generally find the task-agnostic argument compelling, though think this could have been communicated better in the paper. The runtime experiment and discussion is great, but this needs to be in the paper (and the authors have not indicated that they will add it to the paper).
Thank you for your response. I appreciate the additional experiments and the discussion about the proposed framework's limitations.
**W6:** Examples of PRDS p-values include independent p-values and p-values from test statistics that are jointly normally distributed, if all correlations between test statistics are positive (`Wang, Ruodu, and Aaditya Ramdas. "False discovery rate control with e-values." Journal of the Royal Statistical Society Series B: Statistical Methodology 84.3 (2022): 822-852.`).

**W7:** We appreciate these suggestions. We have thoroughly addressed these issues in the revised manuscript.

**Q1:** The “privacy-preserving feature” of PSPS refers to the fact that we only require summary statistics as input for inference, rather than individual-level raw data (features $X$ and label $Y$). This terminology is commonly used in human genetics, healthcare, and multicenter electronic health record analysis to describe methods that enable statistical inference without directly accessing personal data. This approach is analogous to federated learning (`Kairouz, Peter, et al. 2021`) in the machine learning literature. For example, consider a scenario where labeled data is in one center and unlabeled data in another, yet researchers cannot access individual-level data from both centers simultaneously. Under such conditions, current ML-assisted inference, which relies on accessing both labeled and unlabeled data to minimize a joint loss function, is not feasible. However, PSPS circumvents this issue by aggregating summary statistics from multiple centers, thus performing statistical inference while upholding the privacy of individual-level data.

We also acknowledge that in the machine learning literature, "privacy-preserving" often specifically refers to techniques like differential privacy. To avoid confusion, we have revised the terminology from “privacy-preserving” to "federated inference" and provided an example to illustrate this use case more clearly.

**Q2:** First, we apologize for the error that we mistakenly used Figure 2e for Figure 2f due to a typo in the code for making the figure. The correct Figure 2f is attached in the rebuttal pdf. However, this does not change our main results related to statistical efficiency: PSPS is more efficient than PPI++ and other existing methods. As we have previously explained in our comparison of PSPS and PPI++, PSPS leverages information across multiple coordinates in a multi-dimensional parameter estimation task, resulting in asymptotic variances that are less than or equal to those produced by PPI++. Regarding the concern about the choice of tuning, the POP-Inf method, which tunes the element-wise variance, also shows wider confidence intervals than PSPS in linear regression.

**Q3:** Our prediction pipeline comprises two components: prediction for unlabeled data and prediction for labeled data. To predict bone mineral density in unlabeled data, we first selected predictive features by 1) calculating the correlation of bone mineral density with 466 other variables (sample size > 200,000 from UK Biobank) using labeled data and 2) selecting the top 50 variables with the highest correlations as inputs for the SoftImpute algorithm to predict bone mineral density in the unlabeled data.

For the labeled data, we employ a similar approach but incorporate 10-fold cross-validation to prevent overfitting. We select the predictive variables and train the SoftImpute model using 90% of the labeled data. We then perform predictions on the remaining 10% in each fold and repeat this process 10 times across all folds.

We have included these details in the Appendix of the manuscript.

**Final note**: We are excited that you find our method flexible, simple and elegant, and the real data application compelling. If you have any further questions, please do not hesitate to let us know. If our responses have resolved your concerns, we kindly request you to consider increasing your score and championing our paper.
**Q (continued):**
* Even for M-estimation tasks that have already been implemented, PSPS offers the additional advantage of relying solely on summary statistics. The “task-specific derivations” mentioned throughout our paper were not only referring to statistical tasks, but also scientific tasks. Real-world data analysis in any scientific discipline often involves conventions and nuisances that require careful consideration. For example, our work is partly motivated from genome-wide association studies (GWAS). Statistically, GWAS is a linear regression that regresses an outcome (e.g., height) on many genetic variants. While the regression-based statistical foundation is simple, conducting a valid GWAS requires accounting for numerous technical issues, such as sample relatedness (i.e., study participants may be genetically related) and population structure (i.e., unrelated individuals of the same ancestry are both genetically and phenotypically similar, creating confounded associations in GWAS). Sophisticated algorithms and software have been developed to address these complex issues (`Mbatchou et al. 2021`). It will be very challenging if all these important features need to be reimplemented in an ML-assisted GWAS framework. With our PSPS protocol, researchers can utilize existing algorithms and software that are highly optimized for genetic applications to perform ML-assisted GWAS. This adaptability is not just limited to GWAS but is a major feature of our approach across scientific domains. A main result of this paper is that PSPS enables researchers to conduct ML-assisted inference using well-established data analysis pipelines.

**L:**
Following this suggestion, we have conducted experiments to compare the computational efficiency of our method (PSPS) against existing methods. Utilizing a dataset with 500 labeled and 10,000 unlabeled data points, PSPS required 1.62 seconds for linear regression and 8.27 seconds for logistic regression using 200 bootstrap resampling. The computation of one-step de-biasing using summary statistics alone took 0.032 seconds for linear regression and 0.033 seconds for logistic regression. Current methods, which estimate asymptotic variance via the closed form derived by the Central Limit Theorem instead of resampling, ranged from 0.024 to 0.049 seconds for linear regression and 0.032 to 0.077 seconds for logistic regression. Although PSPS is slower due to its resampling nature, the overall runtime remains relatively short.

| Method | Linear regression  | Logistic regression |
| ------------- | ------------- |------------- |
| PSPS | 1.62s | 8.27s |
| PPI  | 0.024s  | 0.032s |
| PPI++  | 0.031s  | 0.077s |
| POP-Inf  | 0.049s |0.034s |

We also note that we designed PSPS to utilize summary statistics, aiming to integrate seamlessly with existing computationally efficient software routinely used in data analysis. For example, Regenie is a software employed in GWAS that allows for fast computation of tens of millions of linear regressions using tens of thousands of samples  (`Mbatchou et al. 2021`). Our protocol involves initially generating summary statistics using such a software tool, followed by their integration, allows high computational efficiency. 

**Final note**: We are excited that you find our work tackling an important problem and appreciate our theoretical study in statistical efficiency. If you have any further questions, please do not hesitate to let us know. If our responses have resolved your concerns, we kindly request you to consider increasing your score and championing our paper.