Source-linked AI summary
A modern maximum-likelihood theory for high-dimensional logistic regression
Pragya Sur, Emmanuel J. Candes
TL;DR
The paper asks whether classical logistic-regression approximations remain reliable when the number of features is a substantial fraction of the sample size. It develops a high-dimensional theory for MLE bias, variability, and LRT behavior, finding that classical inference fails while the new predictions remain accurate in finite samples.
Problem
Classical logistic-regression software continues using approximations whose validity is uncertain when p is not negligible compared with n.
Method
The paper develops a theory for independent-covariate logistic models, using approximate message passing to characterize the MLE and signal strength to adjust inference.
Results
The theory accurately predicts MLE bias, MLE variability, and LRT behavior asymptotically and in finite samples with p = 0.2n.
Takeaways & Limitations
Classical software-based estimates, standard errors, and p-values can provide unreliable inference in higher-dimensional logistic regression.
Takeaways & Limitations
Detailed proofs are not provided in the paper, and extending the results to correlated predictors and other generalized linear models remains future work.
Abstract
from arXiv · showhide
Every student in statistics or data science learns early on that when the sample size largely exceeds the number of variables, fitting a logistic model produces estimates that are approximately unbiased. Every student also learns that there are formulas to predict the variability of these estimates which are used for the purpose of statistical inference; for instance, to produce p-values for testing the significance of regression coefficients. Although these formulas come from large sample asymptotics, we are often told that we are on reasonably safe grounds when $n$ is large in such a way that $n \ge 5p$ or $n \ge 10p$. This paper shows that this is far from the case, and consequently, inferences routinely produced by common software packages are often unreliable. Consider a logistic model with independent features in which $n$ and $p$ become increasingly large in a fixed ratio. Then we show that (1) the MLE is biased, (2) the variability of the MLE is far greater than classically predicted, and (3) the commonly used likelihood-ratio test (LRT) is not distributed as a chi-square. The bias of the MLE is extremely problematic as it yields completely wrong predictions for the probability of a case based on observed values of the covariates. We develop a new theory, which asymptotically predicts (1) the bias of the MLE, (2) the variability of the MLE, and (3) the distribution of the LRT. We empirically also demonstrate that these predictions are extremely accurate in finite samples. Further, an appealing feature is that these novel predictions depend on the unknown sequence of regression coefficients only through a single scalar, the overall strength of the signal. This suggests very concrete procedures to adjust inference; we describe one such procedure learning a single parameter from data and producing accurate inference
1 Introduction
When the number of features is not negligible relative to the sample size, classical logistic-regression approximations fail: MLEs are biased and overly variable, while likelihood-ratio inference is unreliable. The paper develops asymptotically exact predictions for these effects and shows they can be adjusted using signal strength.
- Modern datasets motivate testing whether classical logistic-regression approximations remain valid when p is not negligible compared with n.
- Failures in moderately large dimensions: The MLE systematically overestimates effect magnitudes, with most coordinates satisfying |β̂_j| > |β_j|.Positive effects are overestimated, while negative effects are estimated in the same direction with magnitudes that are too large.
- Failures in moderately large dimensions: For coefficients drawn from N(3, 16), the theory predicts estimates centered around 1.499β_j rather than β_j.
- Failures in moderately large dimensions: The biased MLE can drive estimated disease probabilities toward zero or one, especially as the true probability approaches an endpoint.This creates substantial overestimation or underestimation of the chance of being sick.
- Failures in moderately large dimensions: Classical chi-square likelihood-ratio approximations yield misleading p-values, including nonuniform values with excessive mass near zero.The paper reports that classical results used by software packages are erroneous in higher dimensions and unreliable for inference.
- Our contribution: The proposed theory predicts MLE bias, MLE variability, and LRT distribution asymptotically exactly when n and p grow at a fixed ratio.Its predictions are extremely accurate in finite samples such as p = 0.2n and depend on β through the scalar signal strength γ.
2 Main Results
In the proportional high-dimensional regime, logistic-regression MLEs exhibit bias, inflated variability, and non-chi-square likelihood-ratio statistics, all characterized through dimensionality and signal strength.
- Setting: When p/n → κ and the signal strength γ is fixed, the paper characterizes high-dimensional MLE and LLR behavior for independent Gaussian features.The covariates are scaled as X_i ∼ N(0, n^-1I_p), with lim ||β||²/n = γ².
- When does the MLE exist?: γ > g_MLE(κ) makes the MLE asymptotically nonexistent because a separating hyperplane sends its estimate to infinity.Below the boundary, the probability that the MLE exists converges to one; above it, that probability converges to zero.
- Behavior of the MLE: The MLE coordinates are centered around α⋆β, have asymptotic variance σ⋆², and become decorrelated from the signal after centering.The nonlinear system defining the predictions has a unique solution exactly in the region where the MLE asymptotically exists.
- Behavior of the MLE: σ⋆ increases with both dimensionality κ and signal strength γ, while α⋆ approaches 1 as κ approaches zero.Thus, the high-dimensional variance exceeds the classical low-dimensional benchmark, whereas the classical unbiasedness limit is recovered when p = o(n).
- Null coordinates: Null MLE coordinates are asymptotically normal and asymptotically independent for any finite collection of null variables.Their variance is determined by the nonlinear system and increases as the magnitude of the other regression coefficients increases.
- Likelihood-ratio testing: The LRT converges to a rescaled chi-square rather than χ², with a factor larger than one for κ > 0; the approximation is accurate even in distributional tails.Empirical p-value calculations show close agreement, including essentially perfect coverage at n = 8000.
3 Adjusting Inference by Estimating the Signal Strength
The paper estimates the single unknown signal-strength parameter by probing the MLE existence frontier, then plugs it into the theory to adjust inference and predictions.
- Estimating signal strength: The proposed ProbeFrontier method estimates γ by locating the dimensionality where the MLE fails to exist in half of repeated subsamples.It checks MLE existence through linear-programming feasibility across a fine grid and linearly interpolates the 50% transition.
- Estimating signal strength: ProbeFrontier identifies the transition by finding the smallest grid value with estimated nonexistence probability at least 0.5, then setting γ̂ = g_MLE(κ̂).The smooth phase-transition boundary makes interpolation more precise with a sufficiently fine grid.
- Empirical performance: Estimated signal strength remains extremely close to the diagonal across γ values from 0.3 to 5, supporting reasonably precise estimation.Plugging γ̂ into the nonlinear system yields estimates of the bias, standard deviation, and LRT rescaling factor.
- Adjusted inference: Across 6000 replicates, estimated parameters recover the true values to the first digit, and adjusted p-values agree with the uniform distribution in bulk and extreme tails.The result supports practical use of the estimated multiplicative correction for inference.
- Debiasing predictions: Dividing coefficients and predictions by the estimated bias removes the observed extreme shrinkage and yields predictions centered around the regression function.In the example, α⋆ = 1.499 and ProbeFrontier gives α̂ = 1.511; the debiased predictions no longer falsely indicate almost certain outcomes.
4 Main Mathematical Ideas
The paper analyzes the high-dimensional logistic-regression MLE through an AMP algorithm and rigorous state-evolution characterization. This framework yields asymptotic descriptions of MLE behavior and the LRT when the relevant iterations converge.
- AMP framework: An AMP algorithm tracks the logistic-regression MLE as both n and p become large.The algorithm is initialized with β0 and generates iterates {S_t, β̂_t}.
- State evolution: The variance-map updates converge to a unique fixed point (α⋆, σ⋆, λ⋆) in the regime where the MLE exists.Initializing at this fixed point makes the parameter sequence stationary.
- AMP framework: If the AMP iterates converge to a fixed point, that point satisfies the MLE optimality condition ∇ℓ(β̂⋆) = 0.This connects the algorithm’s fixed points directly to the MLE.
- State evolution: The AMP sequence is analyzed by reducing it to generalized AMP, which characterizes limiting variances through state evolution.The asymptotic description uses a standard Gaussian variable, the coefficient distribution Π, and σ⋆.
- MLE distribution: Theorem 2 characterizes AMP-iterate limits using expectations over Π and an independent standard Gaussian variable.A smoothed approximation handles the discontinuous indicator used to study null coordinates.
- LRT distribution: The null-coordinate distribution can be used to derive the LRT distribution, whose high-dimensional form deviates from χ2 because of inflated MLE variance and λ⋆.The latter reflects the non-classical Hessian behavior and eigenvalue spread.
5 Broader Implications and Future Directions
The paper argues that classical logistic-regression inference can fail badly when dimensionality or signal strength is large. It develops a theory for independent-covariate models and identifies correlated predictors and other generalized linear models as future extensions.
- Implications: Classical theory can seriously overestimate effect magnitudes in high-dimensional logistic regression, making risk predictions problematic.The issue worsens as dimensionality, signal strength, or both increase.
- Implications: Classical p-values can be misleading, so routinely used tools may fail for both estimation and testing.This motivates methods that adjust inference in high-dimensional settings.
- Contributions: The new theory describes the asymptotic distributions of the MLE and LRT for models with independent covariates and supports inference adjustment using estimated signal strength.The authors state that the method works extremely well, while noting interest in studying alternatives.
- Future directions: Corresponding results for correlated predictors and extensions to other generalized linear models remain future work.These are identified as two promising directions.
A Fisher information
This appendix defines the Fisher-information-based variance used by classical logistic-regression inference. It shows that, for a null coefficient with γ2 = 5, the predicted standard deviation is 2.66.
- Definition: The Fisher information matrix is defined for the logistic-regression model with Gaussian covariates Xi ∼ N(0, n^-1I).Its entries are obtained from the model’s second-derivative structure.
- Classical variance: The classical variance of β̂j is obtained from the inverse Fisher information matrix.The appendix derives this quantity through a direct calculation.
- Classical variance: 2.66 is the predicted standard deviation for a null coefficient when γ2 = 5.The corresponding expression is ν^-1/2.
- Software inference: Software packages use an approximate Fisher information evaluated at the MLE, treating covariates as fixed and replacing β with β̂.The resulting standard errors support routine statistical inference.
B Properties of fixed points of the AMP algorithm
The fixed-point analysis connects the AMP algorithm to the logistic-regression MLE. Algebra using the proximal mapping shows that an AMP fixed point satisfies the MLE score equation.
- Fixed-point equations: At a fixed point of the AMP updates, the pair (β̂⋆, S⋆) satisfies the fixed-point equations.These equations relate the fitted linear predictor to the proximal mapping.
- Fixed-point equations: The proximal-mapping identity implies Xβ̂⋆ = prox_{λ⋆ρ}(λ⋆y + S⋆).This algebraic relation is the key step linking the AMP fixed point to the likelihood equations.
- MLE connection: Therefore, the fixed point β̂⋆ obeys ∇ℓ(β̂⋆) = 0, the optimality condition for the MLE.The AMP fixed point is thus an MLE candidate under the stated connection.
C Refined analysis of the distribution of a null coordinate
The refined analysis uses leave-one-out arguments to characterize the joint distribution of fitted values and derive the asymptotic behavior of a null MLE coordinate. This complements AMP, which describes bulk behavior but is insufficient for an individual coordinate with arbitrarily many nulls.
- Motivation: AMP describes the MLE’s bulk behavior, while leave-one-out analysis is needed for a single null coordinate when arbitrarily many nulls are present.The leave-one-out construction exploits independence after removing both one predictor and one observation.
- Score-equation reduction: Subtracting the full and reduced score equations yields equations for the null coefficient and the difference between full- and reduced-model estimates.Solving these equations after separating the null coordinate produces an expression for its MLE.
- Asymptotic reduction: The denominator in the null-coordinate expression is asymptotically κ/λ[−j] + oP(1), linking the coordinate calculation to the reduced-model curvature.Here λ[−j] is defined through the trace of the inverse reduced Hessian.
- Leave-one-out approximation: For a null predictor, the reduced-model fitted values can be approximated because the predictor is independent of the response and reduced-model estimator.Dropping one observation is assumed to have little effect on fitted values when observations are numerous and i.i.d.
- Residual averaging: The residuals become asymptotically independent, allowing averages of squared residuals to converge to their corresponding expectations.The resulting fitted-value approximation uses the proximal mapping relation proxλρ(z) + λρ′(proxλρ(z)) = z.
- Joint distribution: The joint limit is characterized through two quantities whose variances are γ2 and κσ2⋆, with their covariance also entering the calculation.This joint distribution leads to the asymptotic variance of the null MLE coordinate.