SLOE: A Faster Method for Statistical Inference in High-Dimensional Logistic Regression
Year: 2022

Abstract:
Logistic regression remains one of the most widely used tools in applied statistics, machine learning and data science. However, in moderately high-dimensional problems, where the number of features d is a non-negligible fraction of the sample size n, the logistic regression maximum likelihood estimator (MLE), and statistical procedures based the large-sample approximation of its distribution, behave poorly. Recently, Sur and Candès [2019] showed that these issues can be corrected by applying a new approximation of the MLE’s sampling distribution in this highdimensional regime. Unfortunately, these corrections are difficult to implement in practice, because they require an estimate of the signal strength, which is a function of the underlying parameters β of the logistic regression. To address this issue, we propose SLOE, a fast and straightforward approach to estimate the signal strength in logistic regression. The key insight of SLOE is that the Sur and Candès [2019] correction can be reparameterized in terms of the corrupted signal strength, which is only a function of the estimated parameters β̂. We propose an estimator for this quantity and show that dimensionality correction with it is accurate in finite samples. We demonstrate the importance of routine dimensionality correction in the Heart Disease dataset from the UCI repository, and a genomics application using the UK Biobank. Compared to the existing ProbeFrontier heuristic, SLOE is conceptually simpler and orders of magnitude faster, making it suitable for routine use. Additionally, we provide consistency guarantees in the relevant high-dimensional regime.

heading: 1 Introduction
text: Logistic regression is a workhorse in statistics, machine learning, data science, and many applied fields. It is a generalized linear model that models a binary scalar outcome Y ∈ {0, 1} conditional on observed features X ∈ Rd via
E[Y | X = x] = g(β>x), with g(t) := 11+exp(−t) , (1.1)
with the coefficients β fit using observed data. Logistic regression is popular as a scientific tool because the model is often accurate, and comes with well-established statistical inference procedures for quantifying uncertainty about the parameters β and predictions g(β>x) at test inputs x. For example, most statistical software packages not only produce predictions from the model, but also summaries such as confidence intervals (CIs) and p-values that enable practitioners to understand the strength of evidence for the prediction in a quantitative way. These widely adopted estimation and ∗Google Research, Brain Team †Google Health
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
statistical inference routines are based on approximations from large-sample asymptotic theory of the maximum likelihood estimator (MLE) β̂. These approximations come from the limiting distribution of β̂ as the sample size n tends toward infinity, but the number of covariates d remains fixed.
Unfortunately, these standard approximations perform poorly when the number of covariates d is of a similar order to the number of samples n, even when the sample size is large [Fan et al., 2019]. Recently, Sur and Candès [2019] showed that, in this setting, the behavior of β̂ in finite samples is better approximated by its limiting distribution in another regime, where the aspect ratio κ := d/n > 0 is held fixed as both n and d grow to infinity. They show that the estimated coefficients β̂ (and therefore, the predictions) have systematically inflated magnitude and larger variance than the standard approximation predicts. The precise characterization of the limiting distribution of β̂ in Sur and Candès [2019] and Zhao et al. [2020] justify a new approximation, which, in principle, facilitates debiasing estimates and constructing CIs and p-values for parameters and test predictions, alike.
The goal of this paper is to make the application of these groundbreaking results practical. A major barrier to adoption currently is that the calculation of these bias and variance corrections requires knowing the signal strength γ2 := Var(β>X), which is challenging to estimate because it is a function of the unknown parameter vector β. Sur and Candès [2019] proposed a heuristic called ProbeFrontier to estimate this quantity, but this approach is computationally expensive, conceptually complex, and hard to analyze statistically. We propose a simpler estimation procedure. Our approach reparameterizes the problem in terms of the corrupted signal strength parameter η2 := limn→∞Var(β̂
>X) that includes the noise in the estimate β̂. This is more straightforward (though non-trivial) to estimate. We propose the Signal Strength Leave-One-Out Estimator (SLOE), which consistently estimates η2, and show that using this for inference yields more accurate CIs in finite samples. Importantly, SLOE takes orders of magnitude less computation than ProbeFrontier, having similar runtime to the standard logistic regression fitting routine.

heading: 2 Preliminaries
text: In this section, we revisit some fundamentals of statistical inference with logistic regression, review recent advances in high dimensional settings by Sur and Candès [2019] and Zhao et al. [2020], and discuss the implications of this paradigm shift in terms of better characterizing practice.

heading: 2.1 Logistic Regression and Statistical Inference
text: Estimates of β := (β1, · · · , βd) in the logistic regression model are usually obtained through maximum likelihood estimation, by maximizing the empirical log-likelihood
β̂ := argmax β∈Rd
1
n n∑ i=1 Yi log(g(β >Xi)) + (1− Yi) log(1− g(β>Xi)). (2.1)
The log likelihood is concave, and has a unique maximizer whenever the outcomes are not linearly separable in the covariates. We will use logistic regression synonymously with maximum likelihood estimation in the logistic regression model, and call β̂ the MLE.
The large-sample asymptotic statistical theory [Lehmann and Romano, 2005] taught in nearly every university Statistics program characterizes the behavior of the estimated coefficients and predictions in the limit as n → ∞ while holding the number of features d fixed. Under this theory, estimates converge to their true value, β̂
p→ β, and the estimation error β̂ − β and prediction error g(β̂>x)− g(β>x) will be small to observe, unless amplified by a factor of √ n, in which case √ n(β̂ − β) d N(0, I−1β ), where Iβ := E[DβXX>] is the Fisher information matrix,
with Dβ := g(β>X)(1− g(β>X)). Of course, when analyzing real data, data scientists only have access to a finite number of samples, and so this theory serves as an approximation characterization of the behavior expected in practice. If the approximation is good, one can make inferences about the underlying data generating distribution. For example, for δ ∈ (0, 1), we can construct confidence intervals (CIs) that contain the true parameters with probability 1−δ. Propagating the uncertainty to predictions gives CIs for the outcome probabilities that can help contextualize machine learning predictions for users such as scientists or clinicians [Kompa et al., 2021].
The large sample approximation is a good one when the number of predictors is small relative to the number of observations. However, in modern scientific studies and data-driven technology, while many samples may be collected, many features of each sample are also collected. In such settings, the number of features d may be significant compared to the sample size n; then, the large sample characterization can be quite misleading.
We can visualize the extent to which this breakdown in assumptions leads to bias. Figure 1 shows the estimated coefficients and calibration curve for a logistic regression model fit on simulated data with n = 4000 and p = 800. Despite n being reasonably large, these illustrate significant departures from the behavior predicted by the large sample asymptotics. The magnitude of the estimated coefficients is systematically overestimated and they are noisy. Additionally, the model is overconfident in its predictions. Despite this, the classification accuracy of the resulting model on the in-distribution test data is 0.70, which is close to the optimal accuracy of 0.77.
Recently, Sur and Candès [2019], Mai et al. [2019], and Zhao et al. [2020] showed that the behavior of the MLE in high-dimensional settings is better characterized by a different asymptotic approximation. In this approximation, the aspect ratio κ = d/n, instead of the number of predictors d, remains fixed as n grows to infinity. Formally, they study a sequence of problems in which n, d → ∞; d/n→ κ > 0, γ2 = Var(β>X) is fixed as d grows, 3 and the MLE exists asymptotically for (κ, γ). Sur and Candès [2019], Zhao et al. [2020], and our work assume Xi is marginally Gaussian, and are complementary to the mixture-of-Gaussians generative model studied by Mai et al. [2019].
The semantics of this asymptotic regime are a bit more complicated than the classical one; we can no longer think of the argument in terms of a fixed problem with increasing data. However, The characterization from this asymptotic theory better reflects the finite-sample behavior of the MLE β̂ when the data has a moderate aspect ratio κ.
For simplicity, we state the key result from this work in the case where the covariates are drawn from a spherical Gaussian, i.e., X ∼ N(0, Id), but note that these results are generalized in Zhao et al. [2020] to arbitrary covariance structures.4 So far, proof techniques only work for normally distributed features, however empirical results (Fig. 4b and previous work) the work for sub-Gaussian features. In this regime, the distribution of the MLE β̂ converges to a normal distribution asymptotically
3This prevents the expectations g(β>X) from converging to 0 or 1 as n grows. To satisfy this, the effect of each covariate shrinks as n grows. It is complementary to sparsity models, being more realistic for many problems where predictive performance comes from pooling many weak signals.
4Simulations suggest that our methods generalize to models including an intercept term when using the intercept correction conjectured in Zhao et al. [2020].
centered around the inflated parameter αβ, for some α > 1. In particular, for any coefficient βj , √ n ( β̂j − αβj ) d N(0, σ2?), (2.2)
and for the predicted logit of an input5 x ∈ Rd with asymptotically finite norm and √ nx>β = O(1),
√ n‖x‖−12 ( β̂>x− αβ>x ) d N(0, σ2?), (2.3)
for constants α and σ?, to be discussed later, that are only functions of κ and γ. The general convergence result when X has an arbitrary covariance is stated in Theorem 3.1 of Zhao et al. [2020].
If α and σ? were known, this theory provides corrected approximate CIs and p-values. For example, taking Φ−1 to be the quantile function of a standard normal distribution:
(1− δ) CI for βj : ( β̂j α ± Φ−1(1− δ/2) σ? αn1/2 ) (2.4)
p-value under βj = 0 : p = Φ−1(Z) where Z = β̂j/ σ?
αn1/2 . (2.5)
Figure 2 shows how the p-values obtained based on these adjustments have the expected behavior (i.e., they are uniform when the null is true), whereas those from the standard approximation do not.

heading: 2.2 Calculating Dimensionality Corrections
text: Sur and Candès [2019] show that the constants α and σ? are determined by the solution to a system of three non-linear equations, which depend implicitly on the aspect ratio κ and the signal strength γ2 := Var(β>X):  κ2σ2? = E [ 2g(Q1)(λg(proxλG(Q2))) 2 ] , 0 = E [g(Q1)Q1λg(proxλG(Q2))] , 1− κ = E [ 2g(Q1)
1+λg′(proxλG(Q2))
] ,
(2.6)
where proxλG(s) = argmint λG(t) + 1 2 (s− t) 2, for G the anti-derivative of g, and( Q1 Q2 ) ∼ N ( 0, [ γ2 −αγ2 −αγ2 α2γ2 + κσ2? ]) . (2.7)
Here, Q1 is a random variable with the same distribution (asymptotically) as the population logits β>X , and Q2 is a random variable with the same distribution (asymptotically) as the logits with the biased MLE plugged in β̂>X . The auxiliary parameter λ corresponds to the limiting average eigenvalue of the inverse Hessian at the MLE, (1/n ∑n i=1Xig ′(β̂>Xi)X > i ) −1, which is useful for making comparisons to asymptotic approximations from standard theory. Sur and Candès [2019] note that this system has a unique solution so long as (κ, γ) take values at which the MLE exists asymptotically. Once set up, these equations are straightforward to solve with numerical methods.
The key quantity needed to set up these equations is the signal strength parameter γ, which is not directly observed. Thus, to correct for high dimensions in practice requires that γ be estimated. This is difficult, as γ is itself a function of the unknown parameter β, and is the main focus of this paper.
Sur and Candès [2019] suggest a heuristic method, called ProbeFrontier, for estimating γ in practice. The idea is to search for the (sub)sample size n′ < n such that the observations in the subsample are linearly separable. Sub-sampling the data changes the aspect ratio κ without changing the signal strength γ. For a fixed signal strength γ, there is a sharp aspect ratio cutoff κ?(γ) (the “frontier”) above which the data is separable with high probability. Based on these ideas, ProbeFrontier then inverts κ?(γ) to estimate γ from the empirical aspect ratio p/n′. This requires subsampling the data repeatedly at various candidate aspect ratios, and for each subsample, checking whether the data are linearly separable using a linear program. Repeatedly checking the separability is computationally expensive, and the statistical behavior near the frontier makes analysis tricky.
5Note here that x is a sequence in n. Up to a 1/ √ d scale factor, these constraints are satisfied for a fixed
sequences of examples drawn from an iid test dataset with high probability.

heading: 3 A Practical Method for Dimensionality-Corrected Inference
text: As we established above, the key to estimating confidence intervals for β or the predictions g(x>β) is obtain a consistent estimate of γ2 to plug in to Eqs. (2.6) and (2.7) to get estimates of α and σ?. Here, we describe our method for obtaining such an estimate. First, we reparameterize the covariance in (2.7) in terms of a more easily-estimated “corrupted signal strength” parameter, and then, we derive a computationally efficient estimator of it by approximating a consistent leave-one-out estimator.
We implemented this estimator in Python, using scikit-learn to perform the MLE, and numpy / scipy [Harris et al., 2020, Virtanen et al., 2020] for the high dimensional adjustment and inference.

heading: 3.1 Reparameterization
text: We begin by reparameterizing the estimating equations (2.6) and (2.7) in terms of a quantity that is easier to estimate. As it appears in Sur and Candès [2019], (2.7) is written in terms of the signal strength γ2 = limn→∞Var(β>X), but it can also be written in terms of η2 := limn→∞Var(β̂>X). Using the asymptotic properties of β̂ from Proposition 2.1, and Lemmas 2.1, 3.1 from Zhao et al. [2020], the variance has a simple form: η2 = α2γ2 + κσ2?. Thus, we can write (2.7) equivalently as(
Q1 Q2
) ∼ N ( 0, [ η2 − κσstar2 −α ( η2 − κσ2? ) −α ( η2 − κσ2? ) η2 ]) . (3.1)
We can think of η2 as the corrupted signal strength of the predictors, corrupted by the fact that we only have an estimate β̂ of the true signal β in our predictors. Note that Var(β̂>X) = E[β̂>Σβ̂], which suggests a natural strategy for estimate η2 if Σ were known. Shortly, we will define the (corrupted) Signal Strength Leave-One-Out Estimator (SLOE), η̂2SLOE, that estimates η 2 without knowing Σ.
By plugging this estimate from the data along with the aspect ratio κ in to Eqs. (2.6) and (3.1), we can derive estimates α̂, σ̂2?, and λ̂ of the bias α, variance σ 2 ? and hessian inflation λ, respectively. Then, we can construct CIs (2.4) and compute p-values (2.5) using these estimates.

heading: 3.2 SLOE Signal Strength Estimator
text: We now show how to estimate the corrupted signal strength η2 = limn Var(β̂>X) = limn E[β̂>Σβ̂], using a computationally efficient approximation to an ideal, but impractical leave-one-out estimator.
The key challenge in estimating η is that the predictor covariance Σ must also be estimated. When Σ̂ and β̂ are estimated from the same set of observed (Xi)ni=1, these quantities have non-trivial dependence even in large samples, rendering the naïve estimate β̂>Σ̂β̂ asymptotically biased.
This dependence can be mitigated by a leave-one-out (LOO) estimator. In particular, let (β̂−i)ni=1 be the set of n LOO MLEs, each calculated with all of the data except the i-th example. Define the LOO estimator as the empirical variance of the LOO logits (β̂>−iXi) n i=1:
η̂2LOO = 1
n n∑ i=1 (β̂>−iXi) 2 −
( 1
n n∑ i=1 β̂>−iXi
)2 . (3.2)
This estimator eliminates the problematic dependence between β̂ and Σ̂ because β̂−i is independent of Xi. Theorem 1 shows that η̂LOO is consistent for η (see the Supplementary Materials for proof). Theorem 1. Assume that γ and κ are in the range where the MLE exists asymptotically (see Sur and Candès [2019, Theorem 1]), and that for each n, p(n) satisfies limn→∞ p/n = κ. Let Σp be a p × p positive definite matrix with a bounded condition number, for all p, and β be a p-dimensional vector satisfying limn→∞ β>Σpβ = γ2, √ nβ d Π, a distribution with a second
moment and ∑ j β 2 j → EΠ[β2]. Assume that Xi
iid∼ N(0,Σp), and Yi follows the logistic regression model with parameter β given Xi, independent of the other observations. Then, η̂2LOO p→ η2.
Unfortunately, this LOO estimator is impractical, as it requires refitting the model n times. To address this, we define the Signal Strength Leave-One-Out Estimator (SLOE). 6 SLOE replaces each
6In the spirit of the method we leave one ‘O’ out.
LOO logit β̂>−iX with an approximation Si that is a simple update to the full-data logit β >X: H = − ∑ j∈I Xjg ′(β̂>Xj)X > j , Ui = X > i H −1Xi, (3.3)
Si = β̂ >Xi +
Ui
1 + g′(β̂>Xi)Ui (Yi − g(β̂>Xi)),
η̂2SLOE = 1
n n∑ i=1 S2i −
( 1
n n∑ i=1 Si
)2 . (3.4)
To derive this estimator, we use leave-one-out techniques inspired by the theoretical analyses in El Karoui et al. [2013] and Sur and Candès [2019]. Following along the derivation in Sur and Candès [2019] to approximate β̂ − β̂−i, we write out the optimality conditions for β̂ and β̂−i respectively. Let I = {1, . . . , n} be the indices of all the examples and I−i = {1, . . . , i− 1, i+ 1, . . . , n} be the indices of all but the i-th example. Then, by the first order optimality conditions for the MLE,∑
j∈I Xj(Yj − g(β̂>Xj)) = 0, and ∑ j∈I−i Xj(Yj − g(β̂>−iXj)) = 0.
Taking the difference between these two equations gives 0 = Xi(Yi − g(β̂>Xi)) + ∑ j∈I−i Xj(g(β̂ > −iXj)− g(β̂>Xj)). (3.5)
We expect the difference between β̂−i and β̂ to be small, so we can approximate the difference g(β̂>−iXj)− g(β̂>Xj) well with a Taylor expansion of g around β̂>Xj ,
0 ≈ Xi(Yi − g(β̂>Xi)) + ∑ j∈I−i Xjg ′(β̂>Xj)X > j (β̂−i − β̂). (3.6)
When the MLE exists, the Hessian of the log likelihood for the full data H and leave-one-out data, H−i = − ∑ j∈I−i Xjg ′(β̂>Xj)X > j = H +Xig ′(β̂>Xi)Xi,
are invertible. Therefore, we can solve for β̂−i− β̂ to get
β̂−i − β̂ ≈ H−1−iXi(Yi − g(β̂ >Xi)).
Then, we can accurately approximate the term β̂>−iXi as
β̂>Xi +X > i H −1 −iXi(Yi − g(β̂ >Xi)).
To estimate all of the matrix inverses efficiently, we can take advantage of the fact that they are each a rank one update away from the full Hessian. Applying the Sherman-Morrison inverse formula [Sherman and Morrison, 1950] gives
X>i H −1 −iXi =
X>i H −1Xi
1 + g′(β̂>Xi)X>i H −1Xi
.
Therefore, inverting one matrix gives us what we need from all n inverse Hessians. Approximating β̂>−iXi in η̂ 2 LOO with this gives us the SLOE estimator, η̂ 2 SLOE. To show that the SLOE estimator η̂ 2 SLOE is consistent, we show the remainders from the Taylor approximation used to derive SLOE, (3.6), vanish in the limit (see the Supplement for proof). Proposition 2. Under the conditions of Theorem 1, the estimators η̂LOO and η̂SLOE are asymptotically equivalent, η̂LOO = η̂SLOE + oP (1), and therefore, η̂SLOE p→ η.
Remark 1: The majority of the proof of Theorem 1 allows us to derive a rate of convergence, η̂2 − β̂>Σβ̂ = oP (n1/4− ), for any > 0. However, to show that β̂>Σβ̂ → η2, our proof relies on the approximate message passing (AMP) construction of Sur and Candès [2019] and Zhao et al. [2020]. Deriving rates of convergence for these estimators would require significant new technical results on the analysis of AMP, none of which currently quantify convergence rates. In Figure 3, we show that the estimation error in η̂SLOE is dominated by this latter difference in simulation, but the empirical error has the shape of a standard n−1/2 rate. ♦

heading: 4 Simulations
text: 0.1 0.2 0.3 0.4 0.5 0.5
0.6
0.7
0.8
0.9
1.0
Co ve
ra ge
90% CI Coverage
Corrected CI Standard CI
(a) Gaussian features
0.1 0.2 0.3 0.4 0.5 0.5
0.6
0.7
0.8
0.9
1.0
Co ve
ra ge
90% CI Coverage
2 0.1 1.0 5.0
(b) Simulated GWAS features
Figure 4. Coverage of true probabilities with 90% standard and corrected CIs over all test examples and over 1000 simulations, with n = 4000. Curves are terminated when the observed data are separable in more than 10% of simulations and simulations with numerical instabilities near the separable frontier (vertical dotted lines) are dropped.
500 1000 2000 4000 8000 16000 Sample size n
1 10
10 0
10 00
10 00 0 Ti m e (s )
SLOE ProbeFrontier
Figure 5. Comparison of wall clock run time for parameter and variance estimation in ProbeFrontier and SLOE over 100 simulations with κ = 0.2 and γ2 = 1 at various sample sizes.
In this section, we evaluate η̂SLOE in simulated data, where ground truth parameter values are known. We show (1) that the asymptotic properties of η̂SLOE hold in realistic sample sizes; (2) that η̂SLOE can be used to make effective dimensionality corrections even when the covariates X do not have a Gaussian distribution, as assumed in the theory; and (3) that η̂SLOE requires orders of magnitude less computation than the current heuristic ProbeFrontier.
Here, we evaluate approximations in terms of the coverage rate of 90% CIs for the prediction probabilities µ := g(β>x) = P (Y = 1 | X = x). In each simulation, we calculate a 90% CI for µi for each observation i in the corresponding test set, using the standard CI and the corrected CI, and compute the fraction of CIs that contain the true µi. Valid 90% CIs will cover the target quantity 90% of the time. This evaluation checks whether the CIs centering and width are appropriately calibrated, and whether the normal approximation adequately captures the tail behavior of the estimator.
In our simulations, we use a data generating process parameterized by the sample size n, the aspect ratio κ, and signal strength γ2. We draw n examples of covariate vectors X from a given distribution, and, conditional on xi, calculate the µi according to the logistic regression model where β is parameterized by γ as follows,
βj =  2γ/ √ d j ≤ d/8 −2γ/ √ d d/8 < j ≤ d/4
0 otherwise,
and finally, draw Yi ∼ Bernoulli (µi). We consider two different distributions for X . First, we generate X ∼ N(0, Ip), as studied in the theory. Second, following an experiment in Sur and Candès [2019], we draw X from a discrete distribution that models feature distributions encountered in Genome-Wide Association Studies (GWAS). Here, each feature Xij takes values in {0, 1, 2}, according to the Hardy-Weinberg equilibrium with parameter pj varying between covariates in the range [0.25, 0.75]. This setting is not explicitly covered by the theory, but we expect the approximation to work well nonetheless because the features have thin tails.
We show the results of these coverage experiments, with n = 4000, in Figure 4a. In both the Gaussian and simulated GWAS settings, across all practically feasible values of κ and γ, the corrected intervals
computed with η̂SLOE reliably provide 90% coverage. This is in stark contrast to CIs generated under standard theory whose coverage begins to drop off at κ > 0.05 (that is, 20 observations per covariate) for all signal strengths. Notably, in the simulated GWAS setting that is not covered by theory, the corrected intervals remain reliable while the standard CIs perform worse than in the Gaussian setting.
SLOE is accurate in finite samples and computationally efficient. Figure 5 compares the time to perform estimation and inference with both SLOE and ProbeFrontier. For all dataset sizes, SLOE is 10x faster, and for n & 3000, SLOE is 100x faster, which is when speed matters the most.

heading: 5 Applications
text: Here, we present applications of dimensionality corrections with η̂SLOE, emphasizing the real-world impact of the differences between standard and corrected statistical inference in logistic regression. See the Supplementary Materials for an additional experiment on the UCI dataset “German Credit Scoring”.
Prediction of Heart Disease Reliable uncertainty quantification is essential for supporting effective decision-making in clinical applications [Kompa et al., 2021]. We show how the dimensionality corrections substantially changes the uncertainty communicated by CIs for predicting a heart disease diagnosis from physiological and demographic information from the Cleveland Clinic Heart Disease dataset [Detrano et al., 1989]. Kompa et al. [2021] showed that logistic regression predictions are sensitive to bootstrapping this training data, demonstrating that uncertainty estimation is important. 7
The Heart Disease dataset (downloadable from the UCI Machine Learning Repository) has 14 predictors, and 303 observations. Converting categorical variables to dummy variables, taking interactions between sex and pain categorization given its clinical significance [Mehta et al., 2019], and balancing the classes, gives 136 training examples and 20 predictors (κ = 0.15). SLOE estimates a bias inflation factor α̂ of 1.40 for this data, so the logits in the standard MLE will be inflated by a factor of 40%. Correcting these inflated logits significantly changes the center of the CIs for downstream predictions. Figure 6 compares the standard CIs to the corrected CIs for a handful of test examples, along with the true labels. Notice that the standard CIs are overconfident, with the top CI being entirely above 99.8% probability of heart disease, when in fact, there wasn’t. The corrected CIs show that in reality, much higher probabilities of misdiagnosis (i.e., 1 in 58 chance of being wrong, rather than 1 in 140) are consistent with the training data. In clinical settings, the uncorrected CIs could communicate an unwarranted sense of confidence in a diagnosis.
Genomics In medical genomics, it is common to build predictors for clinically relevant binary phenotypes with logistic regression models, using genetic variants as features. Given a set of candidate variants (predictors) identified by genome-wide association studies (GWAS) or other related methods, practitioners often use p-values associated with the coefficients of the variants to interpret these models (e.g. Ding et al., 2020). Here, we present a case study with real data, in which the corrected p-values computed with η̂SLOE maintain their intended interpretation, unlike the standard p-values.
7The nonparametric bootstrap turns out to produce invalid CIs when κ > 0; see the Supplementary Materials.
We study the association of genetic predictors with glaucoma, an eye condition that damages the optic nerve. We build a dataset from UK Biobank [Sudlow et al., 2015] with 127 variants directly associated with glaucoma according to a recent meta-analysis of GWAS [Gharahkhani et al., 2020] and randomly sample 1873 “null” single nucleotide polymorphisms (SNPs) also present in HapMap 3 [International HapMap 3 Consortium et al., 2010], such that each sampled SNP is at least 100K base pairs away from all associated variants and from each other, and that pairwise Pearson correlation in the training set (see below) is no greater than 0.8, to avoid correlated features due to linkage disequilibrium. This generates 2000 genetic predictors in total. Glaucoma status was determined based on touchscreen questionnaire responses and hospital episode statistics as performed previously [Khawaja et al., 2018, Alipanahi et al., 2021]. From this, we construct a training set of 5067 cases with glaucoma and 5067 age- and sex-matched controls (κ = 2000/10134 = 0.197), 8 train a logistic regression model, and compute standard and corrected p-values with it.
Figure 7a shows p-values of the randomly selected variants whose effects we expect to be null. The standard p-values are skewed while the corrected p-values follow the expected uniform distribution closely, as we observed in simulated data (Figure 2). Uniform null p-values are important for controlling the error rate in procedures applied to select important variants, e.g., for resource-intensive follow-up. In Figure 7b, we investigate how these p-values interact with the Benjamini-Hochberg (BH) procedure [Benjamini and Hochberg, 1995] for “discovering” important variants while controlling false discovery rates. To compute an approximate empirical false discovery proportion, we take all variants identified by the BH procedure, and compute the proportion of the null variants (as opposed to known glaucoma-associated variants) in this set. The corrected p-values result in better false discovery proportion calibration across a range of false discovery rate targets (Figure 7b).

heading: 6 Discussion
text: The groundbreaking work by Sur and Candès [2019] has significant implications for theoretical understanding, uncertainty quantification, and generally, practical statistical inference in large scale data. Connecting this theory to practice highlights some interesting insights and questions.
Bias and Variance In the regime that we analyze, these dimensionality corrections constitute a something of a statistical free lunch: the corrections reduce both bias and variance, or, in the case of confidence intervals, they both improve coverage to their nominal level (Figure 4a) while shortening their length on the logit scale (Figure 6b). This occurs because bias is created by overfitting,9 as opposed to underfitting, which we can correct. This systematic bias in the parameters appears in downstream predictions, and unless corrected, does not vanish asymptotically. Because the bias correction shrinks the coefficient estimates by a factor of 1/α, it also produces lower-variance predictions, even after correcting for the too-small variance predicted by standard theory.
CIs and Regularized Predictions CIs with valid coverage are an attractive tool for quantifying uncertainty about predictions, but in practice, point predictions are usually made with models that are regularized to reduce the variance of the predictions, at the cost of some amount of bias. Yet, the CIs
8We reserve 1126 individuals (563 cases, 563 controls) sampled in the same way for a test set. 9To overfit, the model reduces training residuals by amplifying coefficient to produce extreme predictions.
around the predictions that we produce are centered around unbiased estimates of the logits; thus, one might ask whether these can produce coherent summaries together. In simulations, we find that the 1− δ CIs contain the optimal regularized predictor (selected via leave-one-out cross validation) more than a 1− δ fraction of the time. This makes sense, as the LOOCV regularized predictor should be closer to the true probability, which we know is within the CIs 1− δ of the time. If one preferred the CIs to always include the regularized predictions, taking the convex hull of the CIs and regularized predictions would give slightly larger, but more interpretable, CIs with higher-than-nominal coverage.
Broader Implications for ML It would be interesting to further explore the implications of these observations in machine learning. The logistic regression MLE is closely related to the cross entropy / softmax loss used frequently in deep learning. The logistic regression model is directly related to softmax probabilities through the independence of irrelevant alternatives assumption, potentially enabling multi-class extensions of this work. However, the paradigm studied here focuses on the underparameterized regime, whereas many deep learning models are overparameterized. Even in large datasets where this may not be true, the signal is often strong enough that that the data are separable, and so the MLE does not exist. Salehi et al. [2019] characterizes generalization into regimes where the data are separable by adding regularization, but their results rely more on assuming isotropic Gaussian features, whereas this work and Zhao et al. [2020] generalize well theoretically to elliptical Gaussian features and experimentally to sub-Gaussian features.
Acknowledgments and Disclosure of Funding
This research has been conducted using the UK Biobank Resource under Application Number 17643. The authors thank Babak Alipanahi and Farhad Hormozdiari for their help and feedback on the genomics experiments, Arjun Seshadri for his feedback on the theoretical results and advice on linear programming to implement ProbeFrontier efficiently, and D. Sculley for his support and feedback on the research.

Review summaries:
The authors propose to make practical the high dimensional analysis of the logistic regression performed by [Sur and Candes]. Indeed, [Sur and Candes] proved that logistic regression in the classical regime (number of samples going to infinity while the dimensionality is fixed) is biased and must be corrected for modern applications where the number of samples and their size are of the same order of magnitude. However, in [Sur and Candes]'s paper, the asymptotic analysis depends on parameters that are difficult to estimate since they depend on ground truth parameters (regression parameter $\beta$). The authors of the present article propose a fast estimation method.

The paper is based on a work by Sur and Candès (2019) who propose a way to correct the bias of the MLE estimator and correct confidence intervals and p-values for logistic regression. The results in Sur Candes are asymptotic results n \to \infty but assume a fixed aspect ration $p / n = \kappa$ where p is the number of regression coefficients. The results have  been refined in Zhao et al (2020),: in Zhao, the regressors are Gaussian but with an arbitrary covariance matrix whereas the regressors are spherical Gaussian in Sur Candes  (the distribution of the regressor plays a key role to compute the correction factor).   The original algorithm by Sur and Candes require to estimate the signal strength, defined as the $\gamma^2 = Var(\beta^T X)$. This is of course challenging because the parameter is unknown. The authors proposed to use the "corrupted signal strength, defined as $\Var(\hat{\beta}^T X)$ where $\hat{\beta}$ is the MLE (it is claimed that the corrupted signal strength has been introduced in Zhao (2020) but on the arxiv version of this paper I did not find such quantity). The correction factors can be computed by solving an equation to involving the corrupted signal strength. The purpose of this paper is to compute an estimator of the corrupted signal strength.