Method-of-Moments Inference for GLMs and Doubly Robust Functionals under Proportional Asymptotics

abstractIn this paper, we consider the estimation of regression coefficients and signal-to-noise (SNR) ratio in high-dimensional Generalized Linear Models (GLMs), and explore their implications in inferring popular estimands such as average treatment effects in high dimensional observational studies. Under the “proportional asymptotic” regime and Gaussian covariates with known (population) covariance $\mathbf{\Sigma}$, we derive $\sqrt{n}$-Consistent and Asymptotically Normal (CAN) estimators of our targets of inference through a Method-of-Moments type of estimators that bypasses estimation of high dimensional nuisance functions and hyperparameter tuning altogether. Additionally, under non-Gaussian covariates, we demonstrate universality of our results under certain additional assumptions on the regression coefficients and $\mathbf{\Sigma}$. We also demonstrate that knowing $\mathbf{\Sigma}$ can be relaxed in our proposed methodology. Finally, we complement our theoretical results with extensive numerical experiments, in comparisons with competing methods.

bibunit[plainnat] \section{Introduction} Statistical inference in Generalized Linear Models (GLMs) nelder1972generalized, mccullagh1989generalized, although a classical topic in statistics, has witnessed renewed enthusiasm in the modern big data era spurred by both theoretical and computational challenges that arise therein jankova2018semiparametric, cai2023statistical, sur2019modernb, sur2019likelihood, candes2020phase, zhao2022asymptotic. This line of research has in turn found resonance in the challenges encountered in the context of inference in observational studies chernozhukov2018double, athey2018approximate, jiang2025new, yadlowsky2022explaining, celentano2023challenges. Specifically, estimation of quantities like the causal effect of an exposure on an outcome or estimation of population quantities under missing data typically relies on understanding nuisance functions such as outcome regression and propensity scores robins1994estimation, scharfstein1999adjusting. These regression functions are often modeled as suitable GLMs, when one needs to adjust for confounders possibly larger in dimension than the available sample size. There now exists a dedicated and comprehensive methodology to deal with inference in both GLMs or observational studies with high dimensional covariates/confounders focused on ideas based on semiparametric theory zhang2014confidence, javanmard2014confidence, jankova2018semiparametric, athey2018approximate, smucler2019unifying, bradic2019minimax, bradic2019sparsity, tan2020model, dukes2021inference, wang2024debiased, liu2023root, su2023estimated. Indeed, this class of methods, in turn, relies on rates of convergence for consistent estimators of high dimensional GLM parameters negahban2012unified. However, even the mere existence of such a consistent estimator relies on a priori unknown low-dimensional (such as sparsity) assumptions in respective GLMs verzelen2012minimax, collier2017minimax, cai2017confidence, bellec2022biasing. To complement the above framework, recent times have witnessed a parallel focus to deal with cases when the entire GLM parameter vector cannot be estimated consistently, and yet there are potential low dimensional summaries of them that can yield themselves to desirable inferential strategies. As a byproduct, one can possibly provide reliable estimation in observational studies. One specific instance that has become popular is when the GLM parameters grow proportionally to the sample size in dimension and do not satisfy additional low-dimensional assumptions bean2013optimal, el2013asymptotic, donoho2016high, lei2018asymptotics\footnote{In the econometrics literature, similar problems have also been studied in (partially) linear models cattaneo2018inference, cattaneo2019two under the name “many-regressor asymptotics” or in settings with many weak instrumental variables newey2009generalized, mikusheva2022inference under the name “many-instrument/many weak IV asymptotics”.}. To reflect the inherent difficulty of this setup in terms of the information-theoretic impossibility of estimating the GLM parameters consistently verzelen2012minimax, barbier2019optimal, recent research has coined it as the “inconsistency regime” celentano2023challenges, and fundamental ideas have already started to carve the contours of this paradigm. A major theme of research in this regime often pertains to initial progress made under Gaussian covariates with known covariance bellec2023debiasing and subsequent demonstration of a universality principle zhao2022asymptotic, dudeja2023universality, han2023universality. Indeed, the Gaussian assumption is not necessarily a simplifying assumption in the development and analysis of the methods in this literature -- one requires highly involved probabilistic machinery to produce sharp analyses of the derived estimators. Nevertheless, the assumption of Gaussianity, coupled with the knowledge or a sufficiently accurate estimator of the population (variance-)covariance matrix of baseline covariates inject enough structure and a priori information to simplify the process. Such structure and information enable one to bypass complicated estimators and analyses, while still achieving remarkably parallel guarantees for inference in GLMs and observational studies in the proportional asymptotic high-dimensional regime. The primary aim of this article is to take steps in that direction. Specifically, we demonstrate that for GLMs with link function meeting certain conditions (see Assumption (ref) later), it is possible to construct a diffeomorphism between functionals of the GLM parameters and carefully crafted low-degree moments of the data, for which $\sqrt{n}$-Consistent and Asymptotic Normal (CAN) estimators exist. This crucial observation forms the core of our proposed methodology. \subsection{Results Highlight} We summarize the main results of the paper below: \begin{enumerate}[label = (\arabic*)] • We propose moments-based identification strategies (for the precise meaning of identification, we refer readers to Lemma (ref) later in the paper) for statistical functionals with nuisance models parameterized as high-dimensional GLMs with the dimension $p$ proportional to $n$ when the covariance matrix of the covariates are known. This allows the construction of estimators of relevant low-dimensional summaries of high-dimensional GLM parameters such as contrasts and Signal-to-Noise Ratio (SNR). Moreover, our methods being reliant on only a few low-dimensional moments of the data are computationally efficient. • Our moment-based identification and estimation strategies generalize to parallel inferential techniques for popular objects of interest in observational studies such as average treatment effects and mean estimands under missing data -- where the analyses depend on two nuisance functions modeled by high dimensional GLMs. Compared to the literature in this class of problems, we do not require sample-splitting and cross-fitting-based ideas owing to our ability to avoid estimating nuisance functions. • Our estimators completely bypass the estimation of high dimensional nuisance parameters and are CAN when the baseline covariates are Gaussian under some additional regularity conditions. We further demonstrate the universality of the proposal beyond Gaussian designs in terms of rates of convergence. • We also demonstrate that the assumption of knowing the population covariance matrix $\mathbf{\Sigma}$ of the design can be dropped for our proposal when the sample covariance matrix estimator of $\mathbf{\Sigma}$ is invertible and $p < c \cdot n$ for some constant $c$, under Gaussian designs. • We conduct extensive numerical experiments to verify the validity of our proposals in finite sample, as well as comparing them with methods from the emerging recent literature. Readers can access the codes for replicating our numerical results from \href{https://github.com/cxy0714/Method-of-Moments-Inference-for-GLMs}{the accompanied GitHub repository}. \end{enumerate} \subsection{Related Works} Our research draws inspiration from several past and ongoing research that aims to address inference in high-dimensional problems. To present a compact survey and comparison with the most related members of this literature we divide our discussions across three broad themes: inference in GLMs, inference for popular observational studies, and the knowledge of variance-covariance matrix of baseline covariates. In each of these sub-parts, we shall further briefly touch upon both ultra-high dimensional regimes under sparsity and proportional high dimensional regimes without sparsity aspects of the results in literature. \subsubsection{Inference in GLMs} In the last two decades, statistical inference for linear and quadratic forms of high dimensional GLM parameters has attracted significant attention from the statistical research community verzelen2012minimax, zhang2014confidence, javanmard2014confidence, dicker2014variance, verzelen2018adaptive, cai2018accuracy, collier2017minimax, guo2022moderate, battey2023inference, celentano2024correlation. Two complementary tracks of emphasis have emerged in this regard. In the first line of activities, the strategy of inference often draws inspiration from classical semiparametric theory jankova2018semiparametric and requires the consistent estimation of ultra high-dimensional (when the dimension $p$ is much larger than the sample size $n$) GLM parameters -- which need to rely on apriori low-dimensional assumptions, such as sparsity, on GLM parameter vectors collier2017minimax, cai2018accuracy, cai2023statistical. To explore regimes where the existence of consistent estimators of entire GLM parameter vectors are impossible, a complementary theme of inference in GLMs has sprung in the last decade under proportional asymptotics (when the dimension $p$ is proportional to the sample size $n$) sur2019modernb, sur2019likelihood, candes2020phase, zhao2022asymptotic. In this regime, the strategy typically involves a careful debiasing surgery on initial suitable yet inconsistent GLM parameter vectors to yield sophisticated CAN estimators of linear and quadratic forms. Indeed, literature in this second direction is more recent and had initially focused on linear models in terms of (i) characterizing the precise risk behavior of convex regularized procedures -- first for Gaussian covariates (see bayati2011lasso, stojnic2013framework, thrampoulidis2018precise, miolane2021distribution, celentano2023lasso and references therein) and then beyond Gaussian gerbelot2020asymptotic, gerbelot2022asymptotic, li2023spectrum, han2023universality; and (ii) inference of linear and quadratic forms of the parameter vector -- albeit mainly in the regime where the design covariance is known apriori bellec2022biasing, bellec2025observable, bellec2023debiasing, song2024hede. Results for GLMs are more complete in terms of estimation of the whole parameter vector using convex regularized methods. Parallel methods in GLMs for inference on linear and quadratic forms are more recent, quite case-specific (e.g. consider binary regression with logistic and/or probit link), do not always cover the whole proportional regime (i.e. all aspects ratio considerations of $p/n$) without further assumptions, and often yield coverage guarantees in an average sense on individual coordinates of the GLM parameter vector instead of individually across coordinates bellec2025observable\footnote{It is noteworthy that bellec2025observable additionally considered Single-Index Models (SIMs) with an unknown link function. We further discuss possible extensions of our work from GLMs to SIMs in Section (ref).}. Another work related to ours is sawaya2023moment, which also concerns statistical inference for GLM parameters. In particular, under assumptions (1) the link function having certain asymmetry (see Section A.8 of sawaya2023moment for a precise statement) and (2) the covariates $\mathbf{X}$ having zero mean, sawaya2023moment use only moments of $Y$ to estimate certain quantities in the State Evolution system, that characterizes the asymptotic behavior of the maximum likelihood and its convex regularized analog, to conduct inference -- thus obviating the requirement of knowing the population covariance matrix $\mathbf{\Sigma}$ of $\mathbf{X}$ or estimating $\mathbf{\Sigma}$ with sufficiently fast convergence rate. However, this important advantage is at the expense of precluding important GLMs such as the logistic or probit regression. Finally, the theoretical results of sawaya2023moment rely on assuming the existence and suitable boundedness of estimators based on minimizing possibly regularized GLM loss functions as well as the existence of unique positive solutions to relevant state evolution equations -- which needs to be further verified and rested outside the scope of the work sawaya2023moment. Since we bypass the estimation of the entire parameter vector while performing CAN estimation of low-dimensional summaries of them, our results do not rely on such further assumptions. \subsubsection{Inference in Observational Studies} Quantities like average treatment effects and mean parameters in missing data problems have now emerged as quintessential examples of functionals in observational studies where the challenges of high dimensional baseline covariates require careful methodological consideration. Similar to the literature in GLM, two complementary themes have emerged here as well -- one regarding ultra-high-dimensional regimes under sparsity and another regarding proportional asymptotic regimes without sparsity but under known Gaussian covariate designs celentano2023challenges or for specific functionals with $p < n$ yadlowsky2022explaining, jiang2025new. Since the ultra-high-dimensional regime under sparsity has been heavily studied athey2018approximate, smucler2019unifying, bradic2019minimax, bradic2019sparsity, tan2020model, tan2020regularized, dukes2021inference, wang2024debiased, liu2023root, the results therein are somewhat complete in terms of necessary and sufficient conditions for CAN estimation. However, without Gaussian covariates or the assumption that $p<n$, neither systematic methods nor CAN guarantees exist for the above examples. Our methods aim to fill this gap in the literature. Finally, we remark that our proposed estimators involve second-order $U$-statistics, thus also drawing connections to the growing literature on using Higher-Order $U$-statistics in semiparametric problems in observational studies robins2008higher, van2021higher, kennedy2024minimax, bonvini2024doubly, breunig2019simple, breunig2024adaptive. Also see Remark (ref) of Section (ref) for a more in-depth discussion. \subsubsection{Known (Population) Covariance} The majority of our results relies on the knowledge of the variance-covariance matrix of baseline covariates in the study. This known (population) covariance assumption has also been consistently imposed in the literature on the inference of high-dimensional GLMs under proportional asymptotics bellec2023debiasing, bellec2025observable, in particular when $p > n$. Indeed, verzelen2018adaptive demonstrate the impossibility of estimating or conducting statistical inference on certain functionals in high dimensional regression with unknown arbitrary variance-covariance matrix of the covariates when $p \gg n^{1 + c}$ for some $c > 0$. This does not preclude the designing of procedures informed by a priori assumptions on the variance-covariance matrix of the covariates -- a philosophy that has indeed been successfully espoused in the ultra-high-dimensional sparse GLM-based inferences verzelen2018adaptive. Its parallels in the proportional asymptotic regime without sparsity assumptions are quite sporadic, and we are only aware of li2023spectrum, and to some extent takahashi2018statistical, that address this problem for linear forms of the coefficients under the right-rotationally invariant design. Our main results also assume that $\mathbf{\Sigma}$ is known. However, under Gaussian designs, in Section (ref), we establish $\sqrt{n}$-consistency of our proposed estimator when $\mathbf{\Sigma}$ is unknown as long as the sample covariance matrix estimator of $\mathbf{\Sigma}$ is invertible, demonstrating that knowing $\mathbf{\Sigma}$ is not essential for our proposal. Furthermore, upon the completion of the first version of our draft, N. Verzelen brought to our attention kong2018estimating. In that paper, the authors developed an estimator of the quadratic form of the regression coefficients (also known as “learnability” in the theoretical computer science literature) in logistic regression only when $\mathbf{X} \sim \mathrm{N}_{p} (\bm{0}, \mathbf{\Sigma})$ with $p \asymp n$, with and without knowing $\mathbf{\Sigma}$. Their estimator without knowing $\mathbf{\Sigma}$ formally resembles the Higher-Order Influence Function estimators robins2008higher, robins2016technical; we will discuss their similarity and difference further in Remark (ref). Our results cover general GLMs beyond logistic regression without forcing the covariates to have zero mean. We also additionally consider more complex functionals often encountered in observational studies, such as the average treatment effects. \subsection*{Organization} To elaborate on the main thesis of the paper, we divide our discussions into the following subsections. In Sections (ref) (knowing $\mathbf{\Sigma}$) and (ref) (not knowing $\mathbf{\Sigma}$), we present our results on inference in GLMs followed by its applications in observational studies collected in Section (ref). Subsequently, Section (ref) validates the theoretical results via numerical experiments. Our article ends by discussing some open problems in Section (ref). All the proof details are deferred to the Appendix. \subsection*{Notation} We denote $(\mathbf{e}_{j}, j = 1, \cdots, p)$ as the standard bases of ${\mathbb{R}}^{p}$. Given any positive integers $\ell \leq k$, we denote $[k] \coloneqq \{1, \cdots, k\}$ and $[\ell:k] \coloneqq \{\ell, \ell + 1, \cdots, k\}$. ${\mathbb{U}}_{n, m} (\cdot)$ is the $m$-th order $U$-statistic operator: given a function $h: {\mathbb{R}}^{m} \rightarrow {\mathbb{R}}$, \begin{align*} {\mathbb{U}}_{n, m} [h (O_{1}, \cdots, O_{m})] \coloneqq \frac{(n - m)!}{n!} \sum_{1 \leq i_{1} \neq \cdots \neq i_{m} \leq n} h (O_{i_{1}}, \cdots, O_{i_{m}}), \end{align*} where $O_{i} \in {\mathbb{R}}$, for $i \in [n]$. When $m = 1$, ${\mathbb{U}}_{n, 1} [h (O)] \equiv n^{-1} \sum_{i = 1}^{n} h (O_{i})$ then reduces to the empirical mean operator. Given any two vectors $\mathbf{v}, \mathbf{u}$ and any matrix ${\mathbb{A}}$ with matching dimensions, we denote $\langle \mathbf{v}, \mathbf{u} \rangle_{{\mathbb{A}}} \coloneqq \mathbf{v}^{\top} {\mathbb{A}} \mathbf{u}$ the inner product between $\mathbf{v}$ and $\mathbf{u}$ with respect to ${\mathbb{A}}$; when ${\mathbb{A}}$ is non-negative semi-definite (n.n.s.d.), given any vector $\mathbf{v}$, this inner product induces a norm $\Vert \mathbf{v} \Vert_{{\mathbb{A}}} \equiv \mathbf{v}^{\top} {\mathbb{A}} \mathbf{v}$. When ${\mathbb{A}} = \mathbf{I}$, the identity matrix, $\Vert \cdot \Vert_{\mathbf{I}} \equiv \Vert \cdot \Vert$ reduces to the standard $\ell_{2}$-norm of a vector. Given a random vector $\mathbf{X}$, $\Vert \mathbf{X} \Vert_{\psi_{2}}$ denotes its Orlicz $\psi_{2}$-norm. $\lambda_{\min} ({\mathbb{A}})$ and $\lambda_{\max} ({\mathbb{A}})$, respectively, denote the minimum and maximum eigenvalues of ${\mathbb{A}}$ when it is symmetric and n.n.s.d.. To avoid clutter, we also introduce the short-hand notation ${\mathbb{E}}^{m} [\cdot] \equiv \{{\mathbb{E}} [\cdot]\}^{m}$. A general theme throughout this paper is to construct a multi-valued map $\Psi = (\Psi_{1}, \cdots, \Psi_{k}): {\mathcal{D}} \rightarrow {\mathcal{R}}$ from its domain ${\mathcal{D}}$ to its range ${\mathcal{R}}$ using MoM. Given a subset $I \subset [k]$, we let $\Psi_{I} \coloneqq \{\Psi_{j}, j \in I\}$ and $\Psi_{I}^{-1}$ be the inverse map of $\Psi_{I}$ if $\Psi_{I}$ is invertible. Given a $k$-th differentiable function $f$, let $f^{(k)}$ denote its $k$-th derivative. Finally, we denote $\Vert f \Vert_{q}$ and $\Vert f \Vert_{\infty}$ as, respectively, the $L_{q} ({\mathbb{P}})$- and $L_{\infty}$-norms of $f$, for $q \geq 1$. \section{Inference in GLMs} In this section, we illustrate our main idea under the following stylized GLM. Suppose that we observe \begin{equation} \tag{$\mathsf{GLM}$} \begin{split} & (Y_i, \mathbf{X}_i)_{i = 1}^{n} \stackrel{\rm i.i.d.}{\sim}\mathbb{P}_{\bm{\beta}}, with Y_i \in \mathbb{R}, \, \mathbf{X}_{i} \sim (\bm{\mu}, \mathbf{\Sigma}) \end{split} \end{equation} where $\bm{\mu} \in {\mathbb{R}}^{p}$ is the unknown mean vector and $\mathbf{\Sigma} \in {\mathbb{R}}^{p \times p}$ is the \textit{known} n.n.s.d. population covariance matrix, and there exists a (possibly) nonlinear \textit{known} link function $\phi: {\mathbb{R}} \rightarrow {\mathcal{R}} \subseteq {\mathbb{R}}$ such that ${\mathbb{E}} [Y | \mathbf{X} = \mathbf{x}] = \phi (\mathbf{x}^\top \bm{\beta})$ with $\bm{\beta} = (\beta_1, \ldots, \beta_p)^{\top} \in \mathbb{R}^p$. The range ${\mathcal{R}}$ of $\phi$ is problem specific -- e.g. when $Y$ is binary and $\phi$ is the expit/logistic function, then ${\mathcal{R}} \equiv [0, 1]$. In this part, we address the question of $\sqrt{n}$-consistent estimation of $\beta_j$ for any $j = 1, \ldots, p$, the linear form of $\bm{\beta}$ along the direction of $\bm{\mu}$ and the quadratic form of $\bm{\beta}$ with respect to $\mathbf{\Sigma}$: \begin{equation} \lambda_{\beta} \coloneqq \bm{\beta}^{\top} \bm{\mu}, \, \, \, \, \gamma_{\beta}^{2} \coloneqq \|\bm{\beta}\|_{\mathbf{\Sigma}}^2. \end{equation} In particular, the quadratic form has been used often in applications related to heritability estimation in genetics, e.g. guo2019optimal, song2024hede and references therein. Moreover, as we will see while studying problems of estimating functionals of interest in observational studies in Section (ref) with two nuisance functions parameterized by GLMs, our analysis in this section will provide the fundamental building blocks. Therefore, looking forward to the case of simultaneously dealing with two high-dimensional GLMs, for the regression coefficients $\bm{\alpha}$ from a separate GLM, we shall also adopt the same convention by denoting $\lambda_{\alpha} \coloneqq \bm{\alpha}^{\top} \bm{\mu}$ and $\gamma_{\alpha}^{2} \coloneqq \Vert \bm{\alpha} \Vert_{\mathbf{\Sigma}}^{2}$. The bounded conditional fourth moment condition on $Y | \mathbf{X}$ is required when establishing the CAN property of our proposed estimators and is imposed here to simplify the exposition (see Appendix (ref)). To discuss the main results of this section and later parts of the paper, we will work with a set of assumptions that we present and discuss before introducing the main ideas of the proposal. It is worth noting that we index assumptions using a single capital letter to highlight their substantive meanings, the majority of which is summarized in Table (ref), together with where the assumptions are imposed throughout this paper. \begin{table}[htbp] \begin{tabular}{c|cc} \hline Assumption & Meaning & Whereabout \\ \hline $\mathsf{D}$ & Bounds on Design Mean & Covariance & Global \\ $\mathsf{L}$ & Link Function & Global \\ $\mathsf{C}$ & Condition Number $p / n$ & Global \\ $\mathsf{B}$ & Bounds on $\Vert \bm{\beta} \Vert$ & Almost Global \\ $\mathsf{V}$ & Conditional Variance/Moments of Response & Almost Global \\ $\mathsf{G}_{0}$ & Gaussian Design and Knowing $\bm{\mu} = \bm{0}$ & Sections (ref) and (ref) \\ $\mathsf{U}_{0}$ & Universality Conditions and Knowing $\bm{\mu} = \bm{0}$ & Sections (ref) and (ref) \\ $\mathsf{G}$ & Gaussian Design with Unknown $\bm{\mu}$ & Except Sections (ref) and (ref) \\ $\mathsf{U}$ & Universality Conditions with Unknown $\bm{\mu}$ & Except Sections (ref) and (ref) \\ \hline \end{tabular} \caption{A glossary for a part of the assumption indices: “Almost Global” means that certain parts of the Assumption are not imposed in some of the theorems. Here $\bm{\mu} \coloneqq {\mathbb{E}} \mathbf{X}$ is the mean of the Design $\mathbf{X}$.} \end{table} First, we state the following “global” assumptions imposed on $\bm{\mu}$, $\mathbf{\Sigma}$, and the link functions $\phi$. \begin{customas}{$\mathsf{D}$} There exist universal constants $M > 0$ such that \begin{align*} M^{-1} \leq \liminf_{p \rightarrow \infty} \lambda_{\min}(\mathbf{\Sigma}) \leq \limsup_{p \rightarrow \infty} \lambda_{\max}(\mathbf{\Sigma}) \leq M \text{ and } \Vert \bm{\mu} \Vert \leq M. \end{align*} \end{customas} \begin{customas}{$\mathsf{L}$} The link function $\phi: {\mathbb{R}} \rightarrow {\mathcal{I}} \subseteq {\mathbb{R}}$, where $\mathcal{I}$ is a closed or open interval in ${\mathbb{R}}$, assumed to satisfy the following conditions: \begin{enumerate} • $\phi$ is three-times differentiable; the first, second, and third derivatives of the link function, together with the link function itself, are integrable with respect to the law of $\mathbf{X}$ and the integrals are all strictly bounded by some universal constant. There also exists a bounded function $f: {\mathbb{R}} \rightarrow {\mathbb{R}}_+$ with $\lim_{|t| \rightarrow \infty} f(t) = 0$ such that $|\phi^{(\ell)} (t)| \leq e^{t^2 f (t)}$ for almost all $t \in {\mathbb{R}}$, for $\ell =1, 2, 3$. • $\phi$ is strictly monotone and both $\phi (x)$ and $\phi^{(1)} (x)$ converge to the boundaries of their respective ranges (possibly $-\infty$ or $+\infty$) as $|x| \rightarrow \infty$. \end{enumerate} \end{customas} \begin{remark} Assumption (ref) accommodates many GLMs commonly encountered in practice, including the logistic regression, probit regression, Poisson/Negative-Binomial log-linear regression, and etc. The latter part of Assumption (ref) (2) also holds for all the above link functions and it will be needed to show that the map from the moments to the functionals of regression coefficients is a global diffeomorphism; see Appendix (ref). \end{remark} As stated in the Introduction, our focus is on estimating functionals related to GLMs within the framework of the proportional asymptotic regime. Consequently, we also operate under the following assumption between the dimension $p$ and the sample size $n$. \begin{customas}{$\mathsf{C}$} There exists $\delta \in (0, \infty)$ such that $\lim_{n \rightarrow \infty} p / n \rightarrow \delta$. \end{customas} \begin{remark} To shorten the exposition, we always assume Assumptions (ref), (ref), and (ref) without explicitly mentioning them unless stated otherwise. \end{remark} Additionally, we state the following boundedness assumption on $\bm{\beta}$; the second part is imposed to rule out the degenerate case $\bm{\beta} \equiv \bm{0}$. Some further comments on this assumption can be found in Remark (ref) later. \begin{customas}{$\mathsf{B}$}\leavevmode \begin{itemize} • There exists a universal constants $0 < \bar{B} < \infty$ such that $\Vert \bm{\beta} \Vert \leq \bar{B}$; • There exists a universal constants $0 < \ubar{B} \leq \bar{B}$ such that $\Vert \bm{\beta} \Vert \geq \ubar{B}$. \end{itemize} \end{customas} Next, we generally need to impose the first part of the following condition on the conditional second moment $Y | \mathbf{X}$. The second part will be needed when establishing the CAN property of our proposed estimator. \begin{customas}{$\mathsf{V}$}\leavevmode \begin{itemize} • We assume that $\Vert \sigma^{2} \Vert_{2}$ is bounded, where $\sigma^{2} (\cdot) \coloneqq {\mathbb{E}} [Y^{2} | \mathbf{X} = \cdot]$ is the conditional second moment function of $Y$ given $\mathbf{X}$; • Let $\sigma^{k} (\mathbf{x}) \coloneqq {\mathbb{E}} [Y^{k} | \mathbf{X} = \mathbf{x}]$ for $k = 2, 4$. We assume that $\Vert \sigma^{k} \Vert_{2}$ is bounded and $\sigma^{k} (\cdot)$ is a GLM sharing the same regression coefficients $\bm{\beta}$, but with possibly different three-times differentiable link functions belonging to $L_{2} ({\mathbb{P}})$. \end{itemize} \end{customas} \begin{remark} Assumption (ref)(2) is mainly made to ease exposition when we establish the CAN property of our proposed estimator (see e.g. Proposition (ref) and Theorem (ref)). But it holds for many popular GLMs encountered in practice. For example, when $Y | \mathbf{X} \sim \mathrm{Ber} (\phi (\bm{\beta}^{\top} \mathbf{X}))$, ${\mathbb{E}} [Y^{2} | \mathbf{X}] = {\mathbb{E}} [Y^{4} | \mathbf{X}] = \phi (\bm{\beta}^{\top} \mathbf{X})$; when $Y | \mathbf{X} \sim \mathrm{Pois} (\phi (\bm{\beta}^{\top} \mathbf{X}))$, ${\mathbb{E}} [Y^{2} | \mathbf{X}] = \phi (\bm{\beta}^{\top} \mathbf{X}) + \phi^{2} (\bm{\beta}^{\top} \mathbf{X})$ and ${\mathbb{E}} [Y^{4} | \mathbf{X}] = \phi^{4} (\bm{\beta}^{\top} \mathbf{X}) + 6 \phi^{3} (\bm{\beta}^{\top} \mathbf{X}) + 7 \phi^{2} (\bm{\beta}^{\top} \mathbf{X}) + \phi (\bm{\beta}^{\top} \mathbf{X})$. \end{remark} Finally, it is worth noting that the symbols for the link function, the regression coefficients and the conditional second moment functions in the above assumptions shall be interpreted as generic notations, as in the sequel we may specialize to problem-specific symbols. For example, later in Section (ref), we also use $\eta$ for the link function and $\bm{\alpha}$ for the regression coefficients. \subsection{Results for designs that are known to have zero mean} To gather intuition for our method, it is instructive first to consider the following assumptions. \begin{customas}{$\mathsf{G_0}$} $\mathbf{X} \sim \mathrm{N}_p (\bm{0}, \mathbf{\Sigma})$ or equivalently $\mathbf{Z} \sim \mathrm{N}_{p} (\bm{0}, \mathbf{I})$ and $\bm{\mu}$ is known to equal $\bm{0}$. \end{customas} Our method then essentially relies on the following result, a direct consequence of Stein's lemma or Gaussian Integration by Parts. \begin{lemma} Under Model (ref), Assumptions (ref), (ref){\rm(1)} and (ref){\rm(1)}, the following hold: \begin{enumerate} • Given any fixed vector $\bm{\upsilon} \in {\mathbb{R}}^p$ and fixed matrix ${\mathbb{M}} \in {\mathbb{R}}^{p\times p}$, the following system of moment equations holds: \begin{subequations} \begin{align} & {\mathbb{E}} [Y \mathbf{X}^{\top}] {\mathbb{M}} {\mathbb{E}} [\mathbf{X} Y] = \bm{\beta}^{\top} \mathbf{\Sigma} {\mathbb{M}} \mathbf{\Sigma} \bm{\beta} \cdot {\mathbb{E}}^2 [\phi' (\mathbf{X}^{\top} \bm{\beta})], \\ & {\mathbb{E}} [Y \mathbf{X}^{\top}] {\mathbb{M}} \bm{\upsilon} = {\mathbb{E}} [\phi'(\mathbf{X}^{\top} \bm{\beta})] \cdot \bm{\beta}^{\top} \mathbf{\Sigma} {\mathbb{M}} \bm{\upsilon}. \end{align} \end{subequations} Consequently, choosing ${\mathbb{M}} = \mathbf{\Sigma}^{-1}$, we have \begin{subequations} \begin{align} & m_{\mathbf{X} Y, 2} \coloneqq {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} Y] = {\mathbb{E}}^2 [\phi' (\mathbf{X}^{\top} \bm{\beta})] \cdot \gamma_{\beta}^{2} = \mathrm{f}_{1}^{2} (\gamma_{\beta}^{2}) \cdot \gamma_{\beta}^{2}, \\ & m_{\mathbf{X} Y, \bm{\upsilon}} \coloneqq {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} \bm{\upsilon} = {\mathbb{E}} [\phi'(\mathbf{X}^{\top} \bm{\beta})] \cdot \bm{\beta}^{\top} \bm{\upsilon} \equiv \mathrm{f}_1 (\gamma_{\beta}^{2}) \cdot \bm{\beta}^{\top} \bm{\upsilon}, \end{align} \end{subequations} where $\mathrm{f}_{1} (t) \coloneqq {\mathbb{E}} [\phi' (Z)]$ with $Z \sim \mathrm{N} (0, t)$ for $t \geq 0$. Denote the map induced by (ref) as $\Psi_{\mathsf{GLM}_{0}, \beta}: \gamma_{\beta}^{2} \mapsto m_{\mathbf{X} Y, 2}$ and the map induced by (ref) as $\Psi_{\mathsf{GLM}_{0}}: (\gamma_{\beta}^{2}, \bm{\beta}^{\top} \bm{\upsilon}) \mapsto (m_{\mathbf{X} Y, 2}, m_{\mathbf{X} Y, \bm{\upsilon}})$. • Further, $\Psi_{\mathsf{GLM}_{0}, \beta}$ is a diffeomorphism with $\nabla (\Psi_{\mathsf{GLM}_{0}, \beta}^{-1})$ bounded; and the same holds for $\Psi_{\mathsf{GLM}_{0}}$. Consequently, $\gamma_{\beta}^{2}$ and $\bm{\beta}^{\top} \bm{\upsilon}$ are identifiable in the sense that the LHS of (ref) uniquely determines the value of $(\gamma_{\beta}^{2}, \bm{\beta}^{\top} \bm{\upsilon})$. \end{enumerate} \end{lemma} We now unpack Lemma (ref), with its proof deferred to the Appendix. Based on the Gaussian design and Stein's lemma, the moment equations (ref) marry the moments on the LHS with certain nonlinear transformation of $\bm{\beta}$ on the RHS. The most important moment equation here is $\Psi_{\mathsf{GLM}_{0}, \beta}$ induced by (ref), that maps the quadratic form $\gamma_{\beta}^{2}$ to the moment $m_{\mathbf{X} Y, 2}$. Assumption (ref) on the link function and Assumption (ref) on $\gamma_{\beta}^{2}$ together ensures that $\Psi_{\mathsf{GLM}_{0}, \beta}$ is a diffeomorphism, so $\gamma_{\beta}^{2} = \Psi_{\mathsf{GLM}_{0}, \beta}^{-1} (m_{\mathbf{X} Y, 2})$ is identified. It will also be made clear later that $\Psi_{\mathsf{GLM}_{0}, \beta}$ being a diffeomorphism entails that $\sqrt{n}$-consistent and CAN estimators of $\gamma_{\beta}^{2}$ can be constructed. After identifying $\gamma_{\beta}^{2}$, by solving (ref), $\bm{\beta}^{\top} \bm{\upsilon} = m_{\mathbf{X} Y, \bm{\upsilon}} / \mathrm{f}_{1} (\gamma_{\beta}^{2})$ is as well identified from the moments. Taking $\bm{\upsilon} = e_{j}$, the $j$-th standard basis in ${\mathbb{R}}^{p}$, the same strategy identifies $\beta_{j}$, for any $j \in [p]$. The conclusions in Lemma (ref), together with all the other identification results under Gaussian designs in this paper, do not require Assumption (ref). However, the above moment equations critically rely on the Gaussianity of $\mathbf{X}$. It is natural to ask if, similar to a growing body of work studying universality for regression models under proportional asymptotics (see e.g. bayati2015universality, montanari2022universality, montanari2023universality, hu2022universality, dudeja2023spectral, lahiry2023universality and references therein), one could move beyond Gaussian designs and demonstrate the universality of the above identification result. We provide a positive answer to this question following a relaxed identification criterion and shifting the burden of assumption from $\mathbf{X}$ to $\bm{\beta}$. \begin{definition}[$\sqrt{n}$-identifiability] We say that a low-dimensional target parameter $\psi \in {\mathbb{R}}^{k}$, where $k$ is strictly bounded, of the underlying statistical model ${\mathbb{P}}$ (e.g. Model (ref)) is $\sqrt{n}$-identifiable if there exists a (possibly) nonlinear map $\Psi$ from $\psi$ to certain moments defined by ${\mathbb{P}}$ induced by $\psi$, such that if given two different values of the target parameter, $\psi$ and $\psi'$, such that $\Vert \psi - \psi' \Vert \gtrsim n^{- 1 / 2}$, then $\Vert \Psi (\psi) - \Psi (\psi') \Vert \gtrsim n^{- 1 / 2}$, for sufficiently large $n$. \end{definition} \begin{customas}{$\mathsf{U_0}$}\leavevmode \begin{enumerate}[label = (\arabic*)] • $\mathbf{X} = \mathbf{\Sigma}^{1 / 2} \mathbf{Z}$, where $\mathbf{Z} = (Z_{1}, \cdots, Z_{p})^{\top}$ has independent coordinates with zero mean and unit variance, and there exists a universal constant $M > 0$ such that $\Vert Z_{j} \Vert_{\psi_{2}} \leq M$ for $j = 1, \cdots, p$; • $\sqrt{p} \mathbf{\Sigma}^{1 / 2} \bm{\beta} \overset{{\mathcal{W}}_{8}}{\rightarrow} \mathsf{b}$ where $\mathsf{b} \sim \rho$ for some probability measure $\rho$ supported on ${\mathbb{R}}$ and we assume that $\rho$ has bounded first and second moments\footnote{Here, given a random vector $\mathsf{A} \in {\mathbb{R}}^{p}$ and a random variable $\mathsf{a} \in {\mathbb{R}}$, the notation $\mathsf{A} \overset{{\mathcal{W}}_{8}}{\rightarrow} \mathsf{a}$ means that the empirical distribution over the coordinates of $\mathsf{A}$ converges in ${\mathcal{W}}_{8}$-distance (8-Wasserstein distance) to the distribution of $\mathsf{a}$, as $p \rightarrow \infty$.}. \end{enumerate} \end{customas} \begin{remark} When $\mathbf{\Sigma} = \mathbf{I}_{p}$, Assumption (ref)(2) reduces to $\sqrt{p} \bm{\beta} \overset{{\mathcal{W}}_{8}}{\rightarrow} \mathsf{b}$. This type of assumptions are commonly imposed . For general population covariance matrix $\mathbf{\Sigma}$, under $\sqrt{p} \bm{\beta} \overset{{\mathcal{W}}_{8}}{\rightarrow} \mathsf{b}$, the more general assumption will be met under additional assumptions on $\mathbf{\Sigma}$. \end{remark} We then have the following parallel result of Lemma (ref), without assuming that $\mathbf{X}$ is Gaussian. \begin{lemma} Under the same assumptions as in Lemma (ref), except with Assumption (ref) replaced by Assumption (ref), the system of moment equations appeared in Lemma (ref) holds approximately with approximation error $O (p^{-3/4}) = O (n^{- 3 / 4})$ as $p \rightarrow \infty$. Thus $\gamma_{\beta}^{2}$ and $\bm{\beta}^{\top} \bm{\upsilon}$ are $\sqrt{n}$-identifiable for any fixed vector $\bm{\upsilon} \in {\mathbb{R}}^{p}$. \end{lemma} The proof of this lemma can be found in Appendix (ref). Taken together, the above results motivate the following MoM-based estimator of $\gamma_{\beta}^{2}$ and $\bm{\beta}^{\top} \bm{\upsilon}$: \begin{equation} \begin{split} & \widehat{\gamma}_{\beta}^{2} \coloneqq \Psi_{\mathsf{GLM}_{0}, \beta}^{-1} \left(\widehat{m}_{\mathbf{X} Y, 2} \mathbbm{1} (\widehat{m}_{\mathbf{X} Y, 2} \in \mathcal{R}_{\mathsf{GLM}_{0}, \beta}) \right), \quad \widehat{m}_{\mathbf{X} Y, 2} \coloneqq {\mathbb{U}}_{n, 2} [Y_1 \mathbf{X}_1^{\top} \mathbf{\Sigma}^{-1} \mathbf{X}_2 Y_2], \\ & \widehat{\bm{\beta}}^{\top} \bm{\upsilon} \coloneqq \frac{\widehat{m}_{\mathbf{X} Y, \bm{\upsilon}}}{\mathrm{f}_{1} (\widehat{\gamma}_{\beta}^{2})}, \quad \widehat{m}_{\mathbf{X} Y, \bm{\upsilon}} \coloneqq {\mathbb{U}}_{n, 1} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} \bm{\upsilon}, \end{split} \end{equation} where $\mathcal{R}_{\mathsf{GLM}_{0}, \beta}$ denotes the range of $\Psi_{\mathsf{GLM}_{0}, \beta}$. Since $\Psi_{\mathsf{GLM}_{0}, \beta}$ is a diffeomorphism, it is clear from the construction above that to prove $\sqrt{n}$-consistency of $\widehat{\gamma}_{\beta}^{2}$ and $\widehat{\bm{\beta}}^{\top} \bm{\upsilon}$ one needs to verify $\max \{\mathsf{var} (\widehat{m}_{\mathbf{X} Y, 2}), \mathsf{var} (\widehat{m}_{\mathbf{X} Y, \bm{\upsilon}})\} = O (1 / n)$, which is indeed the case (see Appendix (ref)). We next summarize the above reasoning as the following proposition. \begin{proposition} Under the Assumptions of Lemma (ref) or Lemma (ref), the following hold: \begin{align*} & \sqrt{n} (\widehat{\gamma}_{\beta}^{2} - \gamma_{\beta}^{2}) = O_{{\mathbb{P}}} (1), \text{ and for any $j = 1, \ldots, p$, } \sqrt{n} (\widehat{\bm{\beta}}^{\top} \bm{\upsilon} - \bm{\beta}^{\top} \bm{\upsilon}) = O_{{\mathbb{P}}} (1). \end{align*} \end{proposition} In fact, as $n \rightarrow \infty$, we can consider a more precise result and record that the above MoM-based estimators are CAN under the Gaussian design and some additional regularity conditions. \begin{proposition} Under Model (ref), Assumptions (ref), (ref) and (ref), if we further assume that $\Vert \bm{\beta} \Vert_{f (\mathbf{\Sigma})}^{2}$ converges to some nontrivial limit for $f (\mathbf{\Sigma}) = \mathbf{\Sigma}, \mathbf{\Sigma}^2, \mathbf{\Sigma}^3$, we have \begin{align*} \sqrt{n} (\widehat{\beta} - \beta_{j}) \overset{\mathcal{L}}{\rightarrow} \mathrm{N} (0, \nu_{j}^{2}) \end{align*} for some constant $\nu_{j}^{2} > 0$ for $j = 1, \cdots, p$ and \begin{align*} \sqrt{n} (\widehat{\gamma}_{\beta}^{2} - \gamma_{\beta}^{2}) \overset{\mathcal{L}}{\rightarrow} \mathrm{N} (0, \nu^{2}) \end{align*} for some constant $\nu^{2} > 0$. \end{proposition} \begin{remark} It is noteworthy that Assumption (ref)(2) is imposed mainly to ensure that the event $\mathbbm{1} (\widehat{m}_{\mathbf{X} Y, 2} \in \mathcal{R}_{\mathsf{GLM}_{0}, \beta})$ holds with probability converging to 1 such that the constraint does not affect the asymptotic distribution of the $U$-statistic estimator $\widehat{m}_{\mathbf{X} Y, 2}$. We conjecture that it is possible to relax this assumption by a more precise analysis of the asymptotic behavior of $\widehat{m}_{\mathbf{X} Y, 2}$ close to the boundary $\Vert \bm{\beta} \Vert = 0$, which is left for future work. \end{remark} \begin{remark} The CAN property of the proposed MoM-based estimators relies on two separate results: (1) the CLT of first-order and second-order $U$-statistics, followed from the results in bhattacharya1992class (see Appendix (ref) for a complete proof) (2) $\Psi_{\mathsf{GLM}_{0}, \beta}$ is a diffeomorphism, and its inverse map, $\Psi_{\mathsf{GLM}_{0}, \beta}^{-1}$, has bounded derivative so the Delta Method can be applied. \end{remark} \begin{remark} In Proposition (ref) (and similar results related to the CAN property of our proposed estimator in the sequel), we need some extra assumptions on the convergence of inner products such as $\bm{\beta}^{\top} f (\mathbf{\Sigma}) \bm{\beta}$. It is worth mentioning that the assumption imposed in the main text might not be tight. By speculating the derivations in Appendix (ref), one only needs either $\bm{\beta}^{\top} \mathbf{\Sigma} \bm{\beta}$ or $\bm{\beta}^{\top} \mathbf{\Sigma}^{3} \bm{\beta}$ to converge. But to avoid unnecessary technical complications that are irrelevant to the main theme of the paper, we decide not to pursue further in this direction. Also, one can easily find sufficient conditions to establish the convergence of such quantities. As a simple example, when $\mathbf{\Sigma} = \mathbf{I}_{p}$ and $\bm{\beta}$ satisfies Assumption (ref), we immediately have $\bm{\beta}^{\top} f (\mathbf{\Sigma}) \bm{\beta}$ to converge. In addition, we also impose an extra condition on ${\mathbb{E}} [Y^{4} | \mathbf{X}]$ via Assumption (ref). This assumption is to ensure that certain re-scaled fourth moments of the $U$-statistics vanish to zero as $n \rightarrow \infty$, which is required based on the proof strategy that we currently employ (see Lemma (ref) and Proposition (ref) in Appendix (ref)). Finally, we do not explicitly specify the form of the asymptotic variance, which is complicated due to the use of second-order $U$-statistics. Nonetheless, in our previous work (Appendix A.8 of liu2024assumption), consistent variance estimators based on tweaking the nonparametric bootstrap have been developed and can be used to conduct inference; also see Section (ref) for a brief discussion and Appendix (ref) for its finite-sample performance. \end{remark} \begin{remark} At this point, readers might wonder why we consider moments such as $m_{\mathbf{X} Y, 2}$ that involves $\mathbf{X}$. When $\bm{\mu} \equiv \bm{0}$, $Y \in \{0, 1\}$ and $\phi (\cdot) = \mathsf{expit} (\cdot)$, it is obvious that all the moments of $Y$ reduce to ${\mathbb{E}} [Y] \equiv 0.5$. Thus without leveraging information of $\mathbf{X}$, it is generally impossible to identify functionals of $\bm{\beta}$. This is why the moment-based estimator in sawaya2023moment, based only on ${\mathbb{E}} [Y]$, cannot be directly applied to logistic or probit regression. \end{remark} \subsection{Results for designs with unknown and possibly non-zero means} Next, we show that the zero covariate-mean condition is not essential to our moment-based approach, but a larger system of moment equations is required to identify relevant parameters in GLMs. The development closely mirrors that of the previous section when $\bm{\mu}$ is known to equal $\bm{0}$. When $\bm{\mu}$ is unknown, however, because the map between the moments and the functionals of the regression coefficients is from ${\mathbb{R}}^{2}$ to ${\mathbb{R}}^{2}$, it is not as straightforward to show that this map is a global diffeomorphism as in the previous section. As before, we first set the stage under Gaussian designs. \begin{customas}{$\mathsf{G}$} $\mathbf{X} \sim \mathrm{N}_p (\bm{\mu}, \mathbf{\Sigma})$ or equivalently $\mathbf{Z} \sim \mathrm{N}_{p} (\bm{0}, \mathbf{I})$. \end{customas} For short, we let $\bm{\nu} = (\nu_{1}, \cdots, \nu_{p})^{\top} \coloneqq \mathbf{\Sigma}^{-1} \bm{\mu}$. The following lemma then generalizes Lemma (ref) to the case of Gaussian designs with unknown mean $\bm{\mu}$. \begin{lemma} Under Model (ref), Assumptions (ref), (ref){\rm(1)}, and (ref){\rm(1)}: \begin{enumerate} • The following system of moment equations holds: \begin{subequations} \begin{align} & m_{Y} \coloneqq {\mathbb{E}} [Y] = \mathrm{f}_{0} (\textcolor{black}{\lambda_{\beta}}, \textcolor{black}{\gamma_{\beta}^{2}}), \\ & m_{\mathbf{X}, 2} \coloneqq {\mathbb{E}} [\mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X}] = \bm{\mu}^{\top} \mathbf{\Sigma}^{-1} \bm{\mu}, \\ & m_{\mathbf{X} Y, \mathbf{X}} \coloneqq {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X}] = m_{Y} \cdot m_{\mathbf{X}, 2} + \mathrm{f}_{1} (\textcolor{black}{\lambda_{\beta}}, \textcolor{black}{\gamma_{\beta}^{2}}) \cdot \textcolor{black}{\lambda_{\beta}}, \\ & m_{\mathbf{X} Y, 2} \coloneqq {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} Y] = m_{Y}^{2} \cdot m_{\mathbf{X}, 2} + \mathrm{f}_{1}^{2} (\textcolor{black}{\lambda_{\beta}}, \textcolor{black}{\gamma_{\beta}^{2}}) \cdot \textcolor{black}{\gamma_{\beta}^{2}} + 2 \cdot m_{Y} \cdot \mathrm{f}_{1} (\textcolor{black}{\lambda_{\beta}}, \textcolor{black}{\gamma_{\beta}^{2}}) \cdot \textcolor{black}{\lambda_{\beta}}, \\ & m_{\nu_{j}} \coloneqq {\mathbb{E}} [\mathbf{X}]^{\top} \mathbf{\Sigma}^{-1} \mathbf{e}_{j} = \bm{\mu}^{\top} \mathbf{\Sigma}^{-1} \mathbf{e}_{j} = \bm{\nu}^{\top} \mathbf{e}_{j} = \nu_{j}, \\ & m_{\beta_{j}} \coloneqq {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} \mathbf{e}_{j} = \mathrm{f}_{0} (\textcolor{black}{\lambda_{\beta}}, \textcolor{black}{\gamma_{\beta}^{2}}) \cdot \nu_{j} + \mathrm{f}_{1} (\textcolor{black}{\lambda_{\beta}}, \textcolor{black}{\gamma_{\beta}^{2}}) \cdot \beta_{j}. \end{align} \end{subequations} where $\mathrm{f}_{k} (s, t) \coloneqq {\mathbb{E}} [\phi^{(k)} (Z)]$ with $Z \sim \mathrm{N} (s, t)$. Denote the forward map induced by this system as \begin{align*} \Psi_{\mathsf{GLM}} = (\Psi_{\mathsf{GLM}, 1}, \Psi_{\mathsf{GLM}, 2}, \cdots \Psi_{\mathsf{GLM}, 6})^{\top}: (\lambda_{\beta}, \gamma_{\beta}^{2})^{\top} \mapsto (m_{Y}, m_{\mathbf{X}, 2}, \cdots, m_{\beta_{j}})^{\top}. \end{align*} • Further, the first four equations of (ref), denoted as $\Psi_{\mathsf{GLM}, [4]}$, can be reduced to \begin{subequations} \begin{align} & m_{1} \coloneqq m_{Y} = \mathrm{f}_{0} (\lambda_{\beta}, \gamma_{\beta}^{2}), \\ & m_{2} \coloneqq m_{\mathbf{X} Y, 2} + m_{Y}^{2} \cdot m_{\mathbf{X}, 2} - 2 \cdot m_{Y} \cdot m_{\mathbf{X} Y, \mathbf{X}} = \mathrm{f}_{1}^{2} (\lambda_{\beta}, \gamma_{\beta}^{2}) \cdot \gamma_{\beta}^{2}. \end{align} Denote the forward map induced by (ref) as $\Psi_{\mathsf{GLM}, \beta} = (\Psi_{\mathsf{GLM}, \beta, 1}, \Psi_{\mathsf{GLM}, \beta, 2})^{\top}: (\lambda_{\beta}, \gamma_{\beta}^{2})^{\top} \mapsto (m_{1}, m_{2})^{\top}$. Then $\Psi_{\mathsf{GLM}, \beta}$ is a diffeomorphism with $\nabla (\Psi_{\mathsf{GLM}, \beta}^{-1})$ bounded. Consequently, $\lambda_{\beta}, \gamma_{\beta}^{2}$, and $\beta_{j}$ are identifiable. \end{subequations} \end{enumerate} \end{lemma} \begin{remark} The proof of Lemma (ref) can be found in Appendix (ref) and Appendix (ref). As mentioned in the beginning of this section, the system of moment equations in (ref) is from ${\mathbb{R}}^{2}$ to ${\mathbb{R}}^{2}$. As a result, showing that (ref) is a global diffeomorphism is nontrivial, which involves the use of Hadamard global inverse function theorem. For details, see Lemma (ref) and Appendix (ref). \end{remark} To establish universality of Lemma (ref) beyond Gaussian designs, we need to first generalize Assumption (ref) to Assumption (ref) below. \begin{customas}{$\mathsf{U}$}\leavevmode \begin{enumerate}[label = (\arabic*)] • $\mathbf{X} = \mathbf{\Sigma}^{1 / 2} \mathbf{Z} + \bm{\mu}$, where $\mathbf{Z} = (Z_{1}, \cdots, Z_{p})^{\top}$ has independent coordinates with zero mean, unit variance, and $\max_{j=1}^p \|Z_j\|_{\psi_2} \leq M$ for some universal constant $M > 0$; • $\sqrt{p} \mathbf{\Sigma}^{1 / 2} \bm{\beta} \overset{{\mathcal{W}}_{8}}{\rightarrow} \mathsf{b}$ and $\sqrt{p} \mathbf{\Sigma}^{- 1 / 2} \bm{\mu} \overset{{\mathcal{W}}_{8}}{\rightarrow} \mathsf{u}$ where $\mathsf{b} \sim \rho$ and $\mathsf{u} \sim \varrho$ respectively for some probability measures $\rho$ and $\varrho$ supported on ${\mathbb{R}}$ and we assume that both $\rho$ and $\varrho$ have bounded first and second moments. \end{enumerate} \end{customas} Lemma (ref) can then be generalized as follows, the proof of which can be found in Appendix (ref). \begin{lemma} Under the same assumptions as in Lemma (ref), except with Assumption (ref) replaced by Assumption (ref), the systems of moment equations appeared in Lemma (ref) hold approximately with approximation error $O (p^{-3 / 4}) = O (n^{-3 / 4})$ as $n \rightarrow \infty$. \end{lemma} When the distribution of $\mathbf{X}$ is unknown and one passes to universality after using a Gaussian identification strategy, we conjecture that certain delocalization conditions, such as Assumption (ref)\rm{(2)} on the (transformed) regression coefficients and covariate mean vector, are necessary. In fact, universality can fail when one starts from a Gaussian identification strategy and Assumption (ref)\rm{(2)} is violated -- as demonstrated via the numerical experiments related to Figures (ref) and (ref); see Section (ref) for more details. \begin{remark} Compared to the special case of knowing $\bm{\mu} \equiv \bm{0}$ in Section (ref), it requires extra moment equations (ref) to (ref) to identify $\gamma_{\beta}^{2}$, together with the linear form $\lambda_{\beta}$. When it is known that $\bm{\mu} = \bm{0}$, (ref), (ref), and (ref), respectively, reduce to constants $1 / 2$, $0$, and $0$. \end{remark} \begin{remark} When the law of $\mathbf{X}$ is absolutely continuous with density $p$ with respect to the Lebesgue measure is (partially) known but non-Gaussian, one could leverage the following generalized Stein's identity (and its higher-order analogues) to obtain similar moment equations \begin{align*} {\mathbb{E}} [f (\mathbf{X}) s (\mathbf{X})] + {\mathbb{E}} [f' (\mathbf{X})] = 0, \end{align*} where $s (\mathbf{x}) \coloneqq \nabla p (\mathbf{x}) / p (\mathbf{x}): {\mathbb{R}}^{p} \rightarrow {\mathbb{R}}^{p}$ is the Stein's score function and $f: {\mathbb{R}}^{p} \rightarrow {\mathbb{R}}$ is any differentiable function such that both terms in the above identity exist. We do not explore these generalizations in this paper and keep it as potential future directions. \end{remark} Gathering the development thus far, we can construct the following estimator of $(\lambda_{\beta}, \gamma_{\beta}^{2})$ based on $\Psi_{\mathsf{GLM}, \beta}$ and its inverse map $\Psi_{\mathsf{GLM}, \beta}^{-1}$: \begin{equation} \begin{split} (\widehat{\lambda}_{\beta}, \widehat{\gamma}_{\beta}^{2}) \coloneqq \Psi_{\mathsf{GLM}, \beta}^{-1} \left( \widehat{m}_{1} \mathbbm{1} \{\widehat{m}_{1} \in {\mathcal{R}}_{\mathsf{GLM}, \beta, 1}\}, \widehat{m}_{2} \mathbbm{1} \{\widehat{m}_{2} \in {\mathcal{R}}_{\mathsf{GLM}, \beta, 2}\} \right), \end{split} \end{equation} where \begin{equation} \begin{split} \widehat{m}_{1} \coloneqq \widehat{m}_{Y} & \coloneqq {\mathbb{U}}_{n, 1} [Y], \quad \widehat{m}_{2} \coloneqq \widehat{m}_{\mathbf{X} Y, 2} + \widehat{m}_{Y}^{2} \cdot \widehat{m}_{\mathbf{X}, 2} - 2 \cdot \widehat{m}_{Y} \cdot \widehat{m}_{\mathbf{X} Y, \mathbf{X}}, \\ \text{and } \widehat{m}_{\mathbf{X}, 2} \coloneqq {\mathbb{U}}_{n, 2} [\mathbf{X}_{1}^{\top} \mathbf{\Sigma}^{-1} \mathbf{X}_{2}], & \, \widehat{m}_{\mathbf{X} Y, \mathbf{X}} \coloneqq {\mathbb{U}}_{n, 2} [Y_{1} \mathbf{X}_{1}^{\top} \mathbf{\Sigma}^{-1} \mathbf{X}_{2}], \, \widehat{m}_{\mathbf{X} Y, 2} \coloneqq {\mathbb{U}}_{n, 2} [Y_{1} \mathbf{X}_{1}^{\top} \mathbf{\Sigma}^{-1} \mathbf{X}_{2} Y_{2}]. \end{split} \end{equation} With $(\widehat{\lambda}_{\beta}, \widehat{\gamma}_{\beta}^{2})$, one can estimate $\beta_{j}$ by simply solving (ref) to (ref): \begin{align*} \widehat{\beta}_{j} \coloneqq \frac{\widehat{m}_{\beta_{j}} - \mathrm{f}_{0} (\widehat{\lambda}_{\beta}, \widehat{\gamma}_{\beta}^{2}) \cdot \widehat{m}_{\nu_{j}}}{\mathrm{f}_{1} (\widehat{\lambda}_{\beta}, \widehat{\gamma}_{\beta}^{2})} \end{align*} where \begin{align*} \widehat{m}_{\nu_{j}} \coloneqq {\mathbb{U}}_{n, 1} [\mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} \mathbf{e}_{j}, \quad \widehat{m}_{\beta_{j}} \coloneqq {\mathbb{U}}_{n, 1} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} \mathbf{e}_{j}. \end{align*} \begin{theorem} Under the Assumptions of Lemma (ref) or Lemma (ref), the following hold: \begin{align*} \sqrt{n} \left( \widehat{\lambda}_{\beta} - \lambda_{\beta} \right) = O_{{\mathbb{P}}} (1), \sqrt{n} \left( \widehat{\gamma}_{\beta}^{2} - \gamma_{\beta}^{2} \right) = O_{{\mathbb{P}}} (1), \text{ and for any $j = 1, \cdots, p$}, \sqrt{n} (\widehat{\beta}_{j} - \beta_{j}) = O_{{\mathbb{P}}} (1). \end{align*} \end{theorem} The final result in this section generalizes Proposition (ref) for the CAN property of our proposed estimators to the case where $\bm{\mu}$ is unknown and possibly non-zero. \begin{theorem} Under Model (ref), Assumptions (ref), (ref) and (ref), if we further assume that $\langle \mathbf{v}_{1}, \mathbf{v}_{2} \rangle_{f (\mathbf{\Sigma})}$ converges to some nontrivial limit for $f (\mathbf{\Sigma}) = \mathbf{\Sigma}^{-1}, \mathbf{\Sigma}, \mathbf{\Sigma}^2, \mathbf{\Sigma}^3$ for $\mathbf{v}_{1}, \mathbf{v}_{2} \in \{\bm{\mu}, \bm{\beta}\}$, we have \begin{align*} \sqrt{n} (\widehat{\beta}_{j} - \beta_{j}) \overset{\mathcal{L}}{\rightarrow} \mathrm{N} (0, \nu_{j}^{2}) \end{align*} for some constant $\nu_{j}^{2} > 0$ for $j = 1, \cdots, p$ and \begin{align*} \sqrt{n} (\widehat{\gamma}_{\beta}^{2} - \gamma_{\beta}^{2}) \overset{\mathcal{L}}{\rightarrow} \mathrm{N} (0, \nu^{2}) \end{align*} for some constant $\nu^{2} > 0$. \end{theorem} In the proof of the above theorem in Appendix (ref), we unpack the conditions of the theorem further. On a higher level, convergence of inner products $\langle \mathbf{v}_{1}, \mathbf{v}_{2} \rangle_{f (\mathbf{\Sigma})}$ to certain nontrivial limits relates to whether the asymptotic variances of moment estimators based on $U$-statistics, after scaled by $\sqrt{n}$, converge to nontrivial limits. For statistical inference, one could construct standard Wald intervals by estimating the asymptotic variances in Theorem (ref) by the bootstrap method developed in liu2024assumption. The performance of bootstrap variance estimators will be assessed in Appendix (ref). \begin{remark}[On the scalings between $p$ and $n$] Even though we focus on the asymptotic regime where $p$ scales proportionally with $n$, it is worth noting that our proposed estimators can achieve consistency even when $p \gg n$, as long as $p = o (n^{2})$. This is because the variances of our $U$-statistic-based moment estimators are generally of order $\frac{1}{n} \vee \frac{p}{n^{2}}$. Asymptotic normality also holds with a different scaling factor $n / \sqrt{p}$ instead of $n^{1 / 2}$, by applying results from bhattacharya1992class. \end{remark} \begin{remark}[Data-driven Approximate Message Passing Schemes for GLMs] Before proceeding, we take a detour to present one immediate application of our MoM-based method for inference in GLMs. Further applications of our method can be found in later Section (ref). For logistic regression, i.e. Model (ref) with $\phi \equiv \mathsf{expit}$, sur2019modernb (and in a more general form, salehi2019impact) showed the following: for $\widehat{\bm{\beta}}_{\mathsf{MLE}} \equiv (\widehat{\beta}_{\mathsf{MLE}, 1}, \cdots, \widehat{\beta}_{\mathsf{MLE}, p})^{\top}$ the Maximum Likelihood Estimator (MLE) of $\bm{\beta}$, letting $(Z_{0}, Z_{1}, Z_{2}, Z_{3})^{\top} \sim \mathrm{N}_{4} (\bm{0}, \mathbf{I}_{4})$, one has \begin{align*} \frac{1}{p} \sum_{j = 1}^{p} h \left( \widehat{\beta}_{\mathsf{MLE}, j} - \bar{\alpha} \beta_{j}, \beta_{j} \right) \overset{{\mathbb{P}}}{\rightarrow} {\mathbb{E}} [h (\bar{\sigma} Z_{0}, {\mathsf{b}})], \end{align*} where $(\bar{\alpha}, \bar{\sigma}, \bar{\gamma})$ is the solution to the following fixed point equations \begin{equation} \begin{split} \left\{ \begin{array}{rl} \delta^{2} \bar{\sigma}^{2} & = \, 2 {\mathbb{E}} \left[ \phi (- \kappa Z_{1}) \left( \bar{\gamma} \phi \left( \mathsf{prox} [\bar{\gamma} \Phi] \left( \kappa \bar{\alpha} Z_{1} + \sqrt{\delta} \bar{\sigma} Z_{2} \right) \right) \right)^{2} \right], \\ 0 & = \, 2 {\mathbb{E}} \left[ \phi (- \kappa Z_{1}) Z_{1} \bar{\gamma} \phi \left( \mathsf{prox} [\bar{\gamma} \Phi] \left( \kappa \bar{\alpha} Z_{1} + \sqrt{\delta} \bar{\sigma} Z_{2} \right) \right) \right], \\ 1 - \delta & = \, 2 {\mathbb{E}} \left[ \phi (- \kappa Z_{1}) \mathsf{prox} [\bar{\gamma} \Phi]' \left( \kappa \bar{\alpha} Z_{1} + \sqrt{\delta} \bar{\sigma} Z_{2} \right) \right]. \end{array} \right. \end{split} \end{equation} Here $\Phi (\cdot) \coloneqq \int_{-\infty}^{\cdot} \phi (t) {\mathrm d} t$ is the anti-derivative of $\phi$, $\kappa \coloneqq {\mathbb{E}} {\mathsf{b}}^{2} \equiv p^{-1} \sum_{j = 1}^{p} {\mathbb{E}} [\beta_{j}^{2}]$ and $\mathsf{prox}$ denotes the proximal operator with $\mathsf{prox} [h]' (\cdot)$ being the first derivative of $\mathsf{prox} [h] (\cdot)$. Since $\kappa$ is unknown, two different methods have been proposed: \texttt{ProbeFrontier} by sur2019modernb leveraging and \texttt{SLOE} by yadlowsky2021sloe using a reparameterization of (ref) in terms of $p^{-1} \sum_{j = 1}^{p} \widehat{\beta}_{\mathsf{MLE}, j}^{2}$ instead of $\kappa$. The proposed MoM-based method evidently offers another approach to estimating $\kappa$ from data. This is also the problem that sawaya2023moment try to solve. However, as mentioned in the Introduction, the approach in sawaya2023moment only applies to asymmetric link functions because it does not rely on $\mathbf{\Sigma}$. \end{remark} \section{The Case of Unknown Population Covariance} In this section, we demonstrate that knowing $\mathbf{\Sigma}$ is not essential to the proposed MoM method in the previous section. To ease presentation, we focus on the case of knowing $\bm{\mu} \equiv \bm{0}$. We will discuss briefly in Appendix (ref) for completeness how to estimate $\mathbf{\Sigma}$ in the general case where both $\bm{\mu}$ and $\mathbf{\Sigma}$ are unknown. In the \href{https://github.com/cxy0714/Method-of-Moments-Inference-for-GLMs}{accompanying GitHub repository}, we also implement the proposed method for this more general scenario. \subsection{A Gaussian-centric approach when \texorpdfstring{$p < n$}{well}} As mentioned in the Introduction, our parameter identification and estimation strategies rely on the knowledge of $\mathbf{\Sigma}$ under the proportional asymptotic regime. In this section, we consider relaxing this assumption under Gaussian designs. To simplify our argument, we assume that $n$ is even and partition all $n$ samples into two equal-sized parts $I_{1} \coloneqq [n / 2]$ and $I_{2} \coloneqq [(n / 2 + 1):2n]$. We use the second partition $I_{2}$ to estimate $\mathbf{\Sigma}$ by the following rescaled sample Gram matrix \begin{align*} \widetilde{\mathbf{\Sigma}} \coloneqq \frac{1}{\frac{n}{2} - p - 1} \sum_{j \in I_{2}} \mathbf{X}_{j} \mathbf{X}_{j}^{\top}. \end{align*} We conjecture that it is possible to attain $\sqrt{n}$-consistency without using sample splitting to estimate $\mathbf{\Sigma}$ liu2023hoif, but we decide to leave the analysis to a future work. When $\mathbf{\Sigma}$ is unknown, we propose to estimate $\gamma_{\beta}^{2}$ as in (ref), except that $\widehat{m}_{\mathbf{X} Y, 2}$ is replaced by \begin{equation} \widehat{m}_{\mathbf{X} Y, 2} \coloneqq \frac{1}{\frac{n}{2} (\frac{n}{2} - 1)} \sum_{i_1 \neq i_2 \in I_{1}} Y_{i_1} \mathbf{X}_{i_1}^{\top} \widetilde{\mathbf{\Sigma}}^{-1} \mathbf{X}_{i_2} Y_{i_2}. \end{equation} It is worth mentioning that, \textit{only in this Section and Appendix (ref)}, $\widehat{m}_{\mathbf{X} Y, 2}$ takes the above form. It is easy to see that $\widehat{m}_{\mathbf{X} Y, 2}$ is an unbiased estimator of $m_{\mathbf{X} Y, 2}$. Moreover, we have the following parallel results to Proposition (ref), the proof of which can be found in Appendix (ref). \begin{proposition} Under Assumptions of Lemma (ref), when $p + 3 < n / 2$ and $p / n \rightarrow c$ for some fixed $c < 1$, the following hold: \begin{align*} \sqrt{n} (\widehat{\gamma}_{\beta}^{2} - \gamma_{\beta}^{2}) = O_{{\mathbb{P}}} (1), \text{ and for any $j = 1, \ldots, p$, } \sqrt{n} (\widehat{\beta}_j - \beta_j) = O_{{\mathbb{P}}} (1). \end{align*} \end{proposition} The additional assumption $p + 3 < n / 2$ is a result of computing the covariance between any two elements of a random matrix drawn from the Inverse-Wishart distribution; see Lemma (ref) for more details. It is worth noting that Proposition (ref) can also be generalized beyond Gaussian designs by using universality-type results on (inverse) sample Gram matrices, e.g. related advances in derezinski2021sparse and derezinski2022unbiased. \begin{remark} Under linear models, one can in fact adapt the estimators of guo2022moderate for unknown $\mathbf{\Sigma}$ by leveraging the linearity structure; for details, see Appendix (ref). Furthermore, when $\mathbf{\Sigma}$ is unknown, without estimating $\mathbf{\Sigma}$, one can still test the null hypothesis $H_{0}: \bm{\beta} = \bm{0}$ by simply estimating the quantity ${\mathbb{E}} [Y \mathbf{X}^{\top}] {\mathbb{E}} [\mathbf{X} Y] = \bm{\beta}^{\top} \mathbf{\Sigma}^{2} \bm{\beta}$, which equals zero if and only if $H_{0}$ is true, using a $U$-statistic ${\mathbb{U}}_{n, 2} (Y_{1} \mathbf{X}_{1}^{\top} \mathbf{X}_{2} Y_{2})$ without the knowledge of $\mathbf{\Sigma}$. \end{remark} \subsection{An alternative approach that works when \texorpdfstring{$p \geq n$}{ill}} The approach considered in the previous section no longer applies to the case of $p \geq n$ (or more precisely, $p \geq n / 2$ when sample splitting is employed), as the sample Gram matrix is not invertible. Here we discuss an alternative approach using Higher-Order $U$-statistics, first considered in kong2018estimating. The idea is to approximate $m_{\mathbf{X} Y, 2} = {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} Y]$, the moment with reciprocal, by Chebyshev polynomials devore1993constructive: \begin{align*} \widetilde{m}_{\mathbf{X} Y, 2, J} \coloneqq \sum_{l = 0}^{J} c_{l} m_{\mathbf{X} Y, 2}^{(l)}, \text{ where } m_{\mathbf{X} Y, 2}^{(l)} \coloneqq {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{l} {\mathbb{E}} [\mathbf{X} Y], \end{align*} up to order $J$ to be specified later. How to specify the Chebyshev polynomial coefficients $\{c_{l}, l = 0, \cdots, J\}$ will be deferred to Appendix (ref). Note that $m_{\mathbf{X} Y, 2}^{(l)}$ can be unbiasedly estimated by the following $(l + 2)$-th order $U$-statistic: \begin{align*} \widehat{m}_{\mathbf{X} Y, 2}^{(l)} \coloneqq \frac{(n - (l + 2))!}{n!} \sum_{1 \leq i_{1} \neq \cdots \neq i_{l + 2} \leq n} Y_{i_{1}} \mathbf{X}_{i_{1}}^{\top} \left\{ \prod_{s = 3}^{l + 2} \mathbf{X}_{i_{s}} \mathbf{X}_{i_{s}}^{\top} \right\} \mathbf{X}_{i_{2}} Y_{i_{2}} \end{align*} with variance of order $\frac{1}{n} \left( \frac{p}{n} \right)^{l + 1}$. Hence we consider the following estimator of $\gamma_{\beta}^{2}$: \begin{align*} \widehat{m}_{\mathbf{X} Y, 2, J (n)} \coloneqq \sum_{l = 0}^{J (n)} c_{l} \widehat{m}_{\mathbf{X} Y, 2}^{(l)}, \end{align*} where $J (n) \asymp (\log n)^{c}$ for some $c$ strictly less than 1. One could approximate the other involved moments and define $\widehat{\gamma}_{\beta, J (n)}^{2}$ and $\widehat{\beta}_{j, J (n)}$ for $j = 1, \cdots, p$ similarly. This alternative approach based on higher-order $U$-statistics also works in principle when $p < n$, but has inflated variance compared to the proposal in the previous section solely based on 2nd-order $U$-statistics. When $p \geq n$, $\widehat{\gamma}_{\beta, J (n)}^{2}$ and $\widehat{\beta}_{j, J (n)}$ are still consistent for $\gamma_{\beta}^{2}$ and $\beta_{j}$ as shown in the Proposition below. The proof is deferred to Appendix (ref). Here we only record the consistency without providing the concrete convergence rate. We conjecture that the tight convergence rate should be slower than $n^{-1 / 2}$ but a more precise characterization demands a much more careful control of the variance of higher-order $U$-statistics and is thus beyond the scope of this article. \begin{proposition} Under Assumptions of Lemma (ref) or Lemma (ref), the following hold: \begin{align*} & \widehat{\gamma}_{\beta, J (n)}^{2} - \gamma_{\beta}^{2} = o_{{\mathbb{P}}} (1), \text{ and for any $j = 1, \ldots, p$, } \widehat{\beta}_{j, J (n)} - \beta_j = o_{{\mathbb{P}}} (1). \end{align*} \end{proposition} \begin{remark} $\widehat{m}_{\mathbf{X} Y, 2}$ defined in (ref) with the prefactor $(\frac{n}{2} - p - 1) / \frac{n}{2}$ removed is nothing but the empirical second-order influence function of the functional ${\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} Y]$, to our knowledge first developed in liu2017semiparametric. In liu2017semiparametric, however, they did not make distributional assumptions on $\mathbf{X}$ such as Assumption (ref) or even Assumption (ref), resulting in more complicated higher-order $U$-statistics for correcting the bias due to estimating $\mathbf{\Sigma}$ by $\widetilde{\mathbf{\Sigma}}$. In particular, the order can diverge to infinity as $n \rightarrow \infty$. Gaussian designs significantly simplify the diverging-order $U$-statistic estimator to a rescaled second-order $U$-statistic. An interesting question to further explore is which other types of assumptions on $\mathbf{X}$, e.g. right rotationally invariant designs li2023spectrum, can also reduce higher-order $U$-statistics to simpler second-order $U$-statistics. When $p \geq n$, here we follow the higher-order $U$-statistic construction of kong2018estimating to avoid estimating $\mathbf{\Sigma}^{-1}$. kong2018estimating also developed a polynomial-time algorithm for computing higher-order $U$-statistics, but only when the $U$-statistic kernels are symmetric. When $\bm{\mu}$ is unknown, the involved $U$-statistics have asymmetric kernels (see Section (ref) and Appendix (ref)), making the computation a much more challenging problem that we plan to address in a separate paper. We mention in passing that robins2016technical also constructed higher-order $U$-statistics, with order diverging at rate $\sqrt{\log n}$, to estimate quantities such as $m_{\mathbf{X} Y, 2}$. Their estimators, however, first estimate the density function of $\mathbf{X}$ and then estimate $\mathbf{\Sigma}^{-1}$ by numerical integration with respect to the estimated density function. Higher-order $U$-statistics are employed to reduce the bias due to density estimation, rather than approximating $\mathbf{\Sigma}^{-1}$ by polynomials. We leave the study of statistical implications of this subtle difference to a future paper. \end{remark} \section{Inference in Observational Studies} In this section, we appeal to the discussions in previous parts to develop estimators for some popular quantities arising in the context of observational studies. We only present theoretical results when the covariates under study have a Gaussian distribution with known covariance $\mathbf{\Sigma}$. Results with unknown $\mathbf{\Sigma}$ can be derived similarly by following the arguments in Section (ref). The universality of the proposed procedure can be established by analogous arguments in the proof of Lemma (ref), and hence omitted to avoid repetition and simplify exposition. The examples we discuss here are popular members of what is now known as the class of Doubly-Robust Functionals robins2008higher, rotnitzky2021characterization, chernozhukov2022locally where in high dimensional instances of the problem practitioners often aim to model the two high dimensional relevant nuisance parameters of the underlying model through flexible GLMs. Although we consider specific and most popular examples in this class, a potential future research avenue is to extend the results in this section to the entire class of Doubly-Robust Functionals characterized in rotnitzky2021characterization. We divide our discussions into three subsections: estimation of causal effects of binary treatment in linear structural models, estimating quantities under missing data, and estimating generalized covariance measure (or equivalently expected conditional covariance liu2024assumption) that have gained popularity in recent literature on conditional independence testing shah2020hardness, christgau2023nonparametric. Throughout this section, we assume that we have access to $n$ i.i.d. copies of the triples $(\mathbf{X}_{i}, A_{i}, Y_{i})_{i = 1}^{n} \overset{\rm i.i.d.}{\sim} {\mathbb{P}}$ where $\mathbf{X}$ is the covariates or the design matrix, except for Section (ref) with a minor modification. The approaches considered in Section (ref) to estimate $\mathbf{\Sigma}^{-1}$ from data apply to all examples here without any conceptual challenges. Consequently, we only present results assuming that $\mathbf{\Sigma}$ is known. \subsection{Causal Effect of a Binary Treatment under Linear Structural Causal Models} First, we consider the task of estimating the causal effect of a binary treatment $A \in \{0, 1\}$ on an outcome $Y$ in a linear structural model. Our goal is to estimate the parameter $\psi$ in the following data generating model with the outcome regression a linear model and the propensity score a GLM with link $\phi$. \begin{equation} \tag{$\mathsf{CE}$} Y = \psi \cdot A + \bm{\beta}^{\top} \mathbf{X} + \varepsilon, \quad A | \mathbf{X} \sim \mathrm{Ber} \left( \phi (\bm{\alpha}^{\top} \mathbf{X}) \right), \end{equation} where $\varepsilon$ has mean zero and bounded second moment given $\mathbf{X}, A$. For convenience, we let $\bm{\mu}_{1} \coloneqq {\mathbb{E}} [\mathbf{X} A]$, $\lambda_{\alpha, 1} \coloneqq \bm{\alpha}^{\top} \bm{\mu}_{1}$ and $\lambda_{\beta, 1} \coloneqq \bm{\beta}^{\top} \bm{\mu}_{1}$. The following lemma characterizes the moment equations under Model (ref), and the derivation can be found in Appendix (ref). \begin{lemma} Under Model (ref), Assumptions (ref), (ref){\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and (ref){\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}, A$, the following system of moment equations holds: \begin{subequations} \begin{align} & m_{A} \coloneqq {\mathbb{E}} [A], \\ & m_{Y} \coloneqq {\mathbb{E}} [Y] = \psi \cdot m_{A} + \lambda_{\beta}, \\ & m_{AY} \coloneqq {\mathbb{E}} [A Y] = (\psi + \lambda_{\beta}) \cdot m_{A} + \gamma_{\alpha, \beta} \cdot \mathrm{f}_{1} (\lambda_{\alpha}, \gamma_{\alpha}^{2}), \\ & m_{\mathbf{X}, 2} \coloneqq {\mathbb{E}} [\mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X}] = \bm{\mu}^{\top} \mathbf{\Sigma}^{-1} \bm{\mu}, \\ & m_{\mathbf{X} A, \mathbf{X}} \coloneqq {\mathbb{E}} [A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X}] = \bm{\mu}_{1}^{\top} \mathbf{\Sigma}^{-1} \bm{\mu}, \\ & m_{\mathbf{X} A, 2} \coloneqq {\mathbb{E}} [A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} A] = m_{A} \cdot m_{\mathbf{X} A, \mathbf{X}} + \mathrm{f}_{1} (\lambda_{\alpha}, \gamma_{\alpha}^{2}) \cdot \lambda_{\alpha, 1}, \\ & m_{\mathbf{X} A, \mathbf{X} Y} \coloneqq {\mathbb{E}} [A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} Y] = \psi \cdot m_{\mathbf{X} A, 2} + \lambda_{\beta, 1}, \\ & m_{\mathbf{X} A Y, \mathbf{X} A} \coloneqq {\mathbb{E}} [Y A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} A] = \psi \cdot m_{\mathbf{X} A, 2} + m_{A} \cdot \lambda_{\beta, 1} + \lambda_{\beta} \cdot m_{A} \cdot m_{\mathbf{X} A, \mathbf{X}} \\ & + \mathrm{f}_{1} (\lambda_{\alpha}, \gamma_{\alpha}^{2}) \cdot \left( \lambda_{\beta} \cdot \lambda_{\alpha, 1} + \gamma_{\alpha, \beta} \cdot m_{\mathbf{X}, \mathbf{X} A} \right) + \mathrm{f}_{2} (\lambda_{\alpha}, \gamma_{\alpha}^{2}) \cdot \gamma_{\alpha, \beta} \cdot \lambda_{\alpha, 1}, \nonumber \end{align} \end{subequations} where the definition of $\mathrm{f}_{k}$ appears in the statement of Lemma (ref). In addition, denote the forward map induced by the RHS of the above system as \begin{align*} \Psi_{\mathsf{CE}} \coloneqq & \left( \Psi_{\mathsf{CE}, 1}, \cdots, \Psi_{\mathsf{CE}, 8} \right): \left( \psi, \lambda_{\alpha}, \gamma_{\alpha}^{2}, \lambda_{\beta}, \gamma_{\alpha, \beta}, \lambda_{\alpha, 1}, \lambda_{\beta, 1} \right) \\ & \mapsto \left( m_{A}, m_{Y}, m_{A Y}, m_{\mathbf{X}, \mathbf{X} A}, m_{\mathbf{X} A, 2}, m_{\mathbf{X} A, \mathbf{X} Y}, m_{\mathbf{X} A Y, \mathbf{X} A} \right). \end{align*} In particular, $\Psi_{\mathsf{CE}}$ is a diffeomorphism. As a consequence, $\psi$, together with $\left( \lambda_{\alpha}, \gamma_{\alpha}^{2}, \lambda_{\beta}, \gamma_{\alpha, \beta}, \lambda_{\alpha, 1}, \lambda_{\beta, 1} \right)$, is identifiable from the above system of moment equations. \end{lemma} The proof of the above lemma can be found in Appendix (ref). As a direct corollary of Lemma (ref), we can construct $\sqrt{n}$-consistent and CAN estimator of $\psi$ as follows. \begin{align*} \widehat{\psi} = \Psi_{\mathsf{CE}, 1}^{-1} \left( \widehat{m}_{A}, \widehat{m}_{Y}, \widehat{m}_{A Y}, \widehat{m}_{\mathbf{X}, \mathbf{X} A}, \widehat{m}_{\mathbf{X} A, 2}, \widehat{m}_{\mathbf{X} A, \mathbf{X} Y}, \widehat{m}_{\mathbf{X} A Y, \mathbf{X} A} \right) \end{align*} where all the moment estimators are constructed in a similar fashion to those in (ref). \begin{theorem} Under Model (ref), Assumptions (ref), (ref){\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and (ref){\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}, A$, we have \begin{align*} \sqrt{n} (\widehat{\psi} - \psi) = O_{{\mathbb{P}}} (1). \end{align*} Alternatively, under Model (ref) with $\varepsilon$ being independent from $\mathbf{X}, A$, Assumptions (ref), (ref) on both $\bm{\alpha}$ and $\bm{\beta}$, and (ref) on both $A | \mathbf{X}$ and $Y | \mathbf{X}, A$, if we further assume that $\langle \mathbf{v}_{1}, \mathbf{v}_{2} \rangle_{f (\mathbf{\Sigma})}$ converges to some nontrivial limit for $f (\mathbf{\Sigma}) = \mathbf{\Sigma}^{-1}, \mathbf{\Sigma}, \mathbf{\Sigma}^2, \mathbf{\Sigma}^3$ for $\mathbf{v}_{1}, \mathbf{v}_{2} \in \{\bm{\mu}, \bm{\alpha}, \bm{\beta}\}$, we have \begin{align*} \sqrt{n} (\widehat{\psi} - \psi) \overset{{\mathcal{L}}}{\rightarrow} \mathrm{N} (0, \nu^{2}) \end{align*} for some constant $\nu^{2} > 0$. \end{theorem} To our knowledge, the above is the first result for CAN estimation for average treatment effect type quantities under the proportional asymptotic regime where $p > n$ is allowed. In this regard, jiang2025new, yadlowsky2022explaining operate under $p < n$ regime, without knowing $\mathbf{\Sigma}$, owing to their reliance on ordinary least squares type estimator for the outcome regression. Our results close this gap in the literature under the knowledge of $\mathbf{\Sigma}$. The proofs of CAN for all the examples in Section (ref), including Theorem (ref), can be found in Appendix (ref). The CAN property of $\widehat{\psi}$ is a corollary of the CAN property established for the estimators of the linear and quadratic forms in Theorem (ref). \subsection{Mean Estimation with Missing Data under Missing-At-Random (MAR)} In this section, we take $A \in \{0, 1\}$ and only observe $Y$ if $A = 1$. Our goal is to estimate $\psi \coloneqq {\mathbb{E}} [Y]$. For $\psi$ to be identifiable from the observed data, the missing data mechanism is assumed to be Missing-At-Random (MAR). More specifically, we assume that ${\mathbb{P}}$ is defined and parameterized as follows: \begin{equation} \tag{$\mathsf{MAR}$} A | \mathbf{X} = \mathbf{x} \sim \mathrm{Ber} \left( \eta (\bm{\alpha}^{\top} \mathbf{x}) \right) \text{ with } \eta (\bm{\alpha}^{\top} \cdot) \in (\ubar{c}, \bar{c}), \, \, Y = \bm{\beta}^{\top} \mathbf{X} + \varepsilon, \end{equation} where $\varepsilon$ has mean zero and bounded second moment and is independent of $A$ and $\mathbf{X}$ and $\ubar{c}, \bar{c}$ are two universal constants satisfying $0 < \ubar{c} < \bar{c} < 1$. As in celentano2023challenges, we assume that only $\mathbf{\Sigma}$ is known. Under Model (ref), \begin{align*} \psi \equiv \bm{\beta}^{\top} \bm{\mu} \equiv {\mathbb{E}} \left[ \frac{A Y}{\eta (\bm{\alpha}^{\top} \mathbf{X})} \right]. \end{align*} This problem is isomorphic to that of estimating treatment specific mean from observational studies under no unmeasured confounding. Let $\mathrm{g}_{j} (t) \coloneqq {\mathbb{E}} [\eta^{(j)} (Z)]$, where $Z \sim \mathrm{N} (\lambda_{\alpha}, \gamma_{\alpha}^{2})$ with $\lambda_{\alpha} \coloneqq \bm{\alpha}^{\top} \bm{\mu}$ and $\gamma_{\alpha}^{2} \coloneqq \Vert \bm{\alpha} \Vert_{\mathbf{\Sigma}}^{2}$. We also define $\gamma_{\alpha, \beta} \coloneqq \bm{\alpha}^{\top} \mathbf{\Sigma} \bm{\beta} = \langle \bm{\alpha}, \bm{\beta} \rangle_{\mathbf{\Sigma}}$. The following lemma characterizes the moment equations under Model (ref), and the derivation can again be found in Appendix (ref). \begin{lemma} Under Model (ref), Assumptions (ref), (ref){\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and (ref){\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, the following system of moment equations holds: \begin{subequations} \begin{align} & m_{A} \coloneqq {\mathbb{E}} [A], \\ & m_{\mathbf{X}, 2} \coloneqq {\mathbb{E}} [\mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X}] = \bm{\mu}^{\top} \mathbf{\Sigma}^{-1} \bm{\mu}, \\ & m_{\mathbf{X} A, \mathbf{X}} \coloneqq {\mathbb{E}} [A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} \bm{\mu} = m_{A} \cdot m_{\mathbf{X}, 2} + \mathrm{f}_{1} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\lambda_{\alpha}} \\ & m_{\mathbf{X} A, 2} \coloneqq {\mathbb{E}} [A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} A] = m_{A}^{2} \cdot m_{\mathbf{X}, 2} + \mathrm{f}_{1}^{2} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\gamma_{\alpha}^{2}} + 2 \cdot m_{A} \cdot \mathrm{f}_{1} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\lambda_{\alpha}} \\ & m_{A Y} \coloneqq {\mathbb{E}} [A Y] = m_{A} \cdot \textcolor{black}{\psi} + \mathrm{f}_{1} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\gamma_{\alpha, \beta}}, \\ & m_{\mathbf{X} A Y, \mathbf{X}} \coloneqq {\mathbb{E}} [Y A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} \bm{\mu} = (m_{A} + m_{\mathbf{X} A, \mathbf{X}}) \cdot \textcolor{black}{\psi} + \left\{ m_{\mathbf{X}, 2} \cdot \mathrm{f}_{1} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) + \mathrm{f}_{2} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\lambda_{\alpha}} \right\} \cdot \textcolor{black}{\gamma_{\alpha, \beta}}, \end{align} \end{subequations} where the definition of $\mathrm{f}_{k}$ appears in the statement of Lemma (ref). In addition, denote the forward map induced by the RHS of the above system as \begin{align*} \Psi_{\mathsf{MAR}} \coloneqq & (\Psi_{\mathsf{MAR}, 1}, \cdots, \Psi_{\mathsf{MAR}, 6}): (\psi, \lambda_{\alpha}, \gamma_{\alpha}^{2}, \gamma_{\alpha, \beta}) \\ & \mapsto (m_{A}, m_{\mathbf{X}, 2}, m_{\mathbf{X} A, \mathbf{X}}, m_{\mathbf{X} A, 2}, m_{AY}, m_{\mathbf{X} A Y, \mathbf{X}}). \end{align*} In particular, $\Psi_{\mathsf{MAR}}$ is a diffeomorphism. As a consequence, $\psi$, together with $\left( \lambda_{\alpha}, \gamma_{\alpha}^{2}, \gamma_{\alpha, \beta} \right)$, is identifiable from the above system of moment equations. \end{lemma} \begin{remark} The system of moment equations for identifying $\psi$ is not unique. For example, one could also replace $m_{\mathbf{X} A Y, \mathbf{X}}$ by $m_{\mathbf{X} A Y, \mathbf{X} A}$, which leads to the following identity: \begin{equation} \tag{11f'} \begin{split} & \ m_{\mathbf{X} A Y, \mathbf{X} A} \coloneqq {\mathbb{E}} [A Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} A] \\ = & \ \left\{ m_{A}^{2} \cdot m_{\mathbf{X}, 2} + 2 m_{A} \cdot \mathrm{f}_{1} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\lambda_{\alpha}} + m_{A}^{2} + \mathrm{f}_{1}^{2} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\gamma_{\alpha}^{2}} \right\} \cdot \textcolor{black}{\psi} \\ & + \left\{ \begin{array}{c} m_{A} \cdot \left( \mathrm{f}_{1} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot m_{\mathbf{X}, 2} + \mathrm{f}_{2} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\lambda_{\alpha}} \right) \\ + \ \mathrm{f}_{1}^{2} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\lambda_{\alpha}} + \mathrm{f}_{1} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \mathrm{f}_{2} (\textcolor{black}{\lambda_{\alpha}}, \textcolor{black}{\gamma_{\alpha}^{2}}) \cdot \textcolor{black}{\gamma_{\alpha}^{2}} \end{array} \right\} \cdot \textcolor{black}{\gamma_{\alpha, \beta}}. \end{split} \end{equation} Combining (ref) with either (ref) or (ref) identifies $\psi$ by solving the corresponding two linear equations. Which combination leads to more efficient estimators of $\psi$ is left for future work to study. If one is also interested in $\mathsf{var} (Y)$, the system (ref) can be further augmented by adding the moment equation for $m_{\mathbf{X} A Y, 2} \coloneqq {\mathbb{E}} [Y A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} A Y]$. \end{remark} The proof of the above lemma can be found in Appendix (ref). As a direct corollary of Lemma (ref), we can construct $\sqrt{n}$-consistent and CAN estimator of $\psi$ as follows. \begin{align*} \widehat{\psi} = \Psi_{\mathsf{MAR}, 1}^{-1} \left( \widehat{m}_{A}, \widehat{m}_{\mathbf{X}, 2}, \widehat{m}_{\mathbf{X} A, \mathbf{X}}, \widehat{m}_{\mathbf{X} A, 2}, \widehat{m}_{AY}, \widehat{m}_{\mathbf{X} A Y, \mathbf{X}} \right) \end{align*} where all the moment estimators are constructed in a similar fashion to those in (ref). \begin{theorem} Under Model (ref), Assumptions (ref), (ref){\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and (ref){\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, we have \begin{align*} \sqrt{n} (\widehat{\psi} - \psi) = O_{{\mathbb{P}}} (1). \end{align*} Alternatively, under Model (ref), Assumptions (ref), (ref) on both $\bm{\alpha}$ and $\bm{\beta}$, and (ref) on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, if we further assume that $\langle \mathbf{v}_{1}, \mathbf{v}_{2} \rangle_{f (\mathbf{\Sigma})}$ converges to some nontrivial limit for $f (\mathbf{\Sigma}) = \mathbf{\Sigma}^{-1}, \mathbf{\Sigma}, \mathbf{\Sigma}^2, \mathbf{\Sigma}^3$ for $\mathbf{v}_{1}, \mathbf{v}_{2} \in \{\bm{\mu}, \bm{\alpha}, \bm{\beta}\}$, we have \begin{align*} \sqrt{n} (\widehat{\psi} - \psi) \overset{{\mathcal{L}}}{\rightarrow} \mathrm{N} (0, \nu^{2}) \end{align*} for some constant $\nu^{2} > 0$. \end{theorem} For estimating $\psi$ in Model (ref), as we mentioned in the Introduction, celentano2023challenges proposed estimators that resemble debiased Lasso in the proportional asymptotic regime, under the Gaussian design with known population covariance matrix $\mathbf{\Sigma}$. We compare the numerical performance of our estimator $\widehat{\psi}$ and those of celentano2023challenges later in Section (ref). \subsection{Estimation of the Generalized Covariance Measure (GCM)} In this section, we assume that ${\mathbb{P}}$ is defined and parameterized as follows: \begin{equation} \tag{$\mathsf{GCM}$} \begin{split} & {\mathbb{E}} [A | \mathbf{X} = \mathbf{x}] = \eta (\bm{\alpha}^{\top} \mathbf{x}), {\mathbb{E}} [Y | \mathbf{X} = \mathbf{x}] = \phi (\bm{\beta}^{\top} \mathbf{x}), \end{split} \end{equation} where both $A$ and $Y$ have bounded second moments. We also let $\sigma_{A}^{2} (\cdot) \coloneqq \mathsf{var} (A | \mathbf{X} = \cdot)$, $\sigma_{Y}^{2} (\cdot) \coloneqq \mathsf{var} (Y | \mathbf{X} = \cdot)$ and $\sigma_{A, Y} (\cdot) \coloneqq \mathsf{cov} (A, Y | \mathbf{X} = \cdot)$. shah2020hardness proposed to test $H_{0}: Y \protect\mathpalette{\protect\independenT}{\perp} A | \mathbf{X}$ by estimating the parameter ${\mathbb{E}} [(Y - \phi (\bm{\beta}^{\top} \mathbf{X})) (A - \eta (\mathbf{X}^{\top} \bm{\alpha}))]$, often referred to as the expected conditional covariance or Generalized Covariance Measure (GCM). Without loss of generality, we take the target of inference as \begin{align*} \psi \coloneqq {\mathbb{E}} [\eta (\bm{\alpha}^{\top} \mathbf{X}) \phi (\mathbf{X}^{\top} \bm{\beta})]. \end{align*} The following lemma characterizes the moment equations under Model (ref), and the derivation can again be found in Appendix (ref). \begin{lemma} Under Model (ref), Assumptions (ref), (ref){\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and (ref){\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, the following system of moment equations holds: \begin{subequations} \begin{align} & m_{A} \coloneqq {\mathbb{E}} [A], \\ & m_{Y} \coloneqq {\mathbb{E}} [Y], \\ & m_{\mathbf{X}, 2} \coloneqq {\mathbb{E}} [\mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X}] = \bm{\mu}^{\top} \mathbf{\Sigma}^{-1} \bm{\mu}, \\ & m_{\mathbf{X} A, \mathbf{X}} \coloneqq {\mathbb{E}} [A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X}] = m_{A} \cdot m_{\mathbf{X}, 2} + \mathrm{g}_{1} (\lambda_{\alpha}, \gamma_{\alpha}^{2}) \cdot \lambda_{\alpha}, \\ & m_{\mathbf{X} Y, \mathbf{X}} \coloneqq {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X}] = m_{Y} \cdot m_{\mathbf{X}, 2} + \mathrm{f}_{1} (\lambda_{\beta}, \gamma_{\beta}^{2}) \cdot \lambda_{\beta}, \\ & m_{\mathbf{X} A, 2} \coloneqq {\mathbb{E}} [A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} A] = m_{A}^{2} \cdot m_{\mathbf{X}, 2} + \mathrm{g}_{1}^{2} (\lambda_{\alpha}, \gamma_{\alpha}^{2}) \cdot \gamma_{\alpha}^{2} + 2 \cdot m_{A} \cdot \mathrm{g}_{1} (\lambda_{\alpha}, \gamma_{\alpha}^{2}) \cdot \lambda_{\alpha}, \\ & m_{\mathbf{X} Y, 2} \coloneqq {\mathbb{E}} [Y \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} Y] = m_{Y}^{2} \cdot m_{\mathbf{X}, 2} + \mathrm{f}_{1}^{2} (\lambda_{\beta}, \gamma_{\beta}^{2}) \cdot \gamma_{\beta}^{2} + 2 \cdot m_{Y} \cdot \mathrm{f}_{1} (\lambda_{\beta}, \gamma_{\beta}^{2}) \cdot \lambda_{\beta}, \\ & m_{\mathbf{X} A, \mathbf{X} Y} \coloneqq {\mathbb{E}} [A \mathbf{X}^{\top}] \mathbf{\Sigma}^{-1} {\mathbb{E}} [\mathbf{X} Y] = m_{A} \cdot m_{Y} \cdot m_{\mathbf{X}, 2} + m_{A} \cdot \mathrm{f}_{1} (\lambda_{\beta}, \gamma_{\beta}^{2}) \cdot \lambda_{\beta} \\ & + m_{Y} \cdot \mathrm{g}_{1} (\lambda_{\alpha}, \gamma_{\alpha}^{2}) \cdot \lambda_{\alpha} + \mathrm{g}_{1} (\lambda_{\alpha}, \gamma_{\alpha}^{2}) \cdot \mathrm{f}_{1} (\lambda_{\beta}, \gamma_{\beta}^{2}) \cdot \gamma_{\alpha, \beta}, \nonumber \end{align} \end{subequations} where the definition of $\mathrm{f}_{k}$ appears in the statement of Lemma (ref). In addition, denote the forward map induced by the RHS of the above system as \begin{align*} \Psi_{\mathsf{GCM}} \coloneqq & (\Psi_{\mathsf{GCM}, 1}, \cdots, \Psi_{\mathsf{GCM}, 8}): (\lambda_{\alpha}, \gamma_{\alpha}^{2}, \lambda_{\beta}, \gamma_{\beta}^{2}, \gamma_{\alpha, \beta}) \\ & \mapsto (m_{A}, m_{Y}, m_{\mathbf{X}, 2}, m_{\mathbf{X} A, \mathbf{X}}, m_{\mathbf{X} Y, \mathbf{X}}, m_{\mathbf{X} A, 2}, m_{\mathbf{X} Y, 2}, m_{\mathbf{X} A, \mathbf{X} Y}). \end{align*} In particular, $\Psi_{\mathsf{GCM}}$ is a diffeomorphism. As a consequence, $\psi$, together with $(\lambda_{\alpha}, \gamma_{\alpha}^{2}, \lambda_{\beta}, \gamma_{\beta}^{2}, \gamma_{\alpha, \beta})$, is identifiable from the above system of moment equations. \end{lemma} The proof of the above lemma can be found in Appendix (ref). As a direct corollary of Lemma (ref), we can construct $\sqrt{n}$-consistent and CAN estimator of $\psi$ as follows. Let \begin{align*} \widehat{\psi} \coloneqq {\mathbb{E}} [\eta (Z_{1}) \phi (Z_{2}], & \text{ where } \left( \begin{matrix} Z_{1} \\ Z_{2} \end{matrix} \right) \sim \mathrm{N}_{2} \left( \left( \begin{matrix} \widehat{\lambda}_{\alpha} \\ \widehat{\lambda}_{\beta} \end{matrix} \right), \left( \begin{matrix} \widehat{\gamma}_{\alpha}^{2} & \widehat{\gamma}_{\alpha, \beta} \\ \widehat{\gamma}_{\alpha, \beta} & \widehat{\gamma}_{\beta}^{2} \end{matrix} \right) \right), \text{ and} \\ (\widehat{\lambda}_{\alpha}, \widehat{\gamma}_{\alpha}^{2}, \widehat{\lambda}_{\beta}, \widehat{\gamma}_{\beta}^{2}, \widehat{\gamma}_{\alpha, \beta}) & \coloneqq \Psi_{\mathsf{GCM}}^{-1} \left( \widehat{m}_{A}, \widehat{m}_{Y}, \widehat{m}_{\mathbf{X}, 2}, \widehat{m}_{\mathbf{X} A, \mathbf{X}}, \widehat{m}_{\mathbf{X} Y, \mathbf{X}}, \widehat{m}_{\mathbf{X} A, 2}, \widehat{m}_{\mathbf{X} Y, 2}, \widehat{m}_{\mathbf{X} A, \mathbf{X} Y} \right), \end{align*} where all the moment estimators are constructed in a similar fashion to those in (ref). To establish $\sqrt{n}$-consistency and CAN of our proposed estimator, we also need to impose the following condition on the conditional covariance function of $A, Y$ given $\mathbf{X}$. \begin{customas}{$\mathsf{Cov}$}\leavevmode \begin{itemize} • We assume that $\Vert \sigma_{A, Y} (\cdot) \Vert_{2}$ is bounded; • We assume that the conditional covariance of $A, Y$ given $\mathbf{X}$ is also a GLM with link function $\zeta$: \begin{align*} \sigma_{A, Y} (\mathbf{X}) = \zeta (\bm{\theta}^{\top} \mathbf{X}) \text{ almost surely,} \end{align*} such that $\zeta$ is three-times differentiable and the first to third derivatives of $\zeta$, together with $\zeta$ itself, are integrable with respect to the law of $\mathbf{X}$ and the integrals are all strictly bounded by some universal constant. \end{itemize} \end{customas} \begin{theorem} Under Model (ref), Assumptions (ref), (ref){\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, (ref){\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, and (ref)\rm{(1)}, we have \begin{align*} \sqrt{n} (\widehat{\psi} - \psi) = O_{{\mathbb{P}}} (1). \end{align*} Alternatively, under Model (ref), Assumptions (ref), (ref) on both $\bm{\alpha}$ and $\bm{\beta}$, (ref) on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, and (ref), if we further assume that $\langle \mathbf{v}_{1}, \mathbf{v}_{2} \rangle_{f (\mathbf{\Sigma})}$ converges to some nontrivial limit for $f (\mathbf{\Sigma}) = \mathbf{\Sigma}^{-1}, \mathbf{\Sigma}, \mathbf{\Sigma}^2, \mathbf{\Sigma}^3$ for $\mathbf{v}_{1}, \mathbf{v}_{2} \in \{\bm{\mu}, \bm{\alpha}, \bm{\beta}, \bm{\theta}\}$, we have \begin{align*} \sqrt{n} (\widehat{\psi} - \psi) \overset{{\mathcal{L}}}{\rightarrow} \mathrm{N} (0, \nu^{2}) \end{align*} for some constant $\nu^{2} > 0$. \end{theorem} Since we do not assume $Y \protect\mathpalette{\protect\independenT}{\perp} A | \mathbf{X}$ and GLMs are only specified for $A | \mathbf{X}$ and $Y | \mathbf{X}$, Assumption (ref)(1) is needed for $\sqrt{n}$-consistency because the variances of the moment estimators based on $U$-statistics involve terms such as ${\mathbb{E}} [A Y \mathbf{X}]$. This is in contrast to the previous two examples. In Section (ref), GLMs are specified for variational independent components $A | \mathbf{X}$ and $Y | \mathbf{X}, A$ of the joint observed-data distribution, whereas in Section (ref), we assume that $Y \protect\mathpalette{\protect\independenT}{\perp} A | \mathbf{X}$. \section{Numerical Experiments} In this section, we verify our theory via extensive numerical experiments. We consider two types of problems. The first problem, studied in Section (ref), is on estimating linear and quadratic forms of regression coefficients $\bm{\alpha}$ of a GLM between the response $A$ and baseline covariates $\mathbf{X}$. The second problem, studied in Section (ref), is on estimating the mean of a response $Y$ subject to missingness under Model (ref). As mentioned, the latter problem is also isomorphic to estimating treatment specific mean from observational studies under no unmeasured confounding. All the results in the numerical experiments are based on 500 Monte Carlos. The dimension-sample size ratio $p / n$ is fixed at 1.2 in Section (ref) and 1.25 in Section (ref). Appendix (ref) contains additional simulation results that are complementary to those reported here in the main text. We always assume that $\mathbf{\Sigma}$ is known except that in Appendix (ref) we conduct simulations to evaluate the performance of our proposed estimators in Section (ref) when knowing neither $\bm{\mu}$ nor $\mathbf{\Sigma}$. The bootstrap variance estimators mentioned right after Theorem (ref) are assessed in Appendix (ref). For all the histograms and normal quantile-quantile plots reported in this section and in Appendix (ref), we choose the numerical experiment with $n = 5000$. The R codes for replicating the numerical experiments can be found in \href{https://github.com/cxy0714/Method-of-Moments-Inference-for-GLMs}{this GitHub repository}. \subsection{Linear and quadratic forms of GLMs} We consider several different experiment settings. But across different settings, the common Data Generating Process (DGP) can be described as follows: $A | \mathbf{X} \sim \mathrm{Ber} \left( \mathsf{expit} \left( \bm{\alpha}^{\top} \mathbf{X} \right) \right)$ where the goal is to estimate a single coordinate $\alpha_{1}$ and the quadratic form $\gamma_{\alpha}^{2} = \bm{\alpha}^{\top} \mathbf{\Sigma} \bm{\alpha}$. Throughout this section, we set the true value of the quadratic form as ${\mathbb{E}}(\gamma_{\alpha}^{2}) \equiv 1$. The true value of $\alpha_{1}$ may vary across different settings. The DGPs vary in the following aspects: \begin{itemize} • $\mathbf{X}_{i} \overset{\rm i.i.d.}{\sim} \mathrm{N}_{p} (\bm{\mu}, \mathbf{\Sigma})$ with $\bm{\mu} \equiv \bm{0}$ and $\mathbf{\Sigma} \equiv \mathbf{I}_{p}$ for $i = 1, \cdots, n$ (Settings 1 & 2); $\mathbf{X}_{i, j} \overset{\rm i.i.d.}{\sim} \mathrm{Rad} (1 / 2)$ for $i = 1, \cdots, n$ and $j = 1, \cdots, p$ where $\mathrm{Rad}$ denotes the Rademacher distribution (Settings 3 & 4). The two cases, though, share the same population mean and covariance matrix. • $\bm{\alpha}$ is dense in the sense that we draw $\bm{\alpha} = (\alpha_{1}, \cdots, \alpha_{p}) \overset{\rm i.i.d.}{\sim} \mathrm{Uniform} ([-\sqrt{3 / p}, \sqrt{3 / p}])$ (Settings 1 & 3); $\bm{\alpha}$ is sparse in the sense that $s = \sqrt{p}$ out of $p$ coordinates are non-zero and each non-zero coefficient equals $p^{- 1 / 4}$ (Settings 2 & 4). Due to the similarity between results for dense and sparse configurations, we defer the figures for the sparse configuration case to Appendix (ref). • We also vary the sample sizes as $n = 1000, 2000, \cdots, 5000$. We note that Setting 2 corresponds to the simulation setting described in Section 6.3 of bellec2025observable. \end{itemize} In each setting, for the problem of estimating certain single coordinate, we compare our MoM-based estimators against the methods of bellec2025observable, which debiases the initial estimator $\widehat{\bm{\alpha}}_{\rm init}$ of $\bm{\alpha}$ using their main Theorem 4.1. We do not further compare different estimators of $\lambda_{\alpha}$ and $\gamma_{\alpha}^{2}$ as the methodology in bellec2025observable is for linear form $\mathbf{v}^{\top} \bm{\alpha}$ with known direction $\mathbf{v}$, not directly applicable here. For simplicity, we use ridge regression to compute $\widehat{\bm{\alpha}}_{\rm init}$. For the tuning parameter $\lambda$ in ridge regression, due to large running time, We choose twelve different values of $\lambda$ ranging from 0.05 to 10, equally spaced on a logarithmic scale. We summarize the results of the numerical experiments below. For single coordinates of $\bm{\alpha}$, we arbitrarily pick $\alpha_{1}$ and $\alpha_{100}$ to present the simulation results. \begin{itemize} • When the design is Gaussian (Settings 1 & 2, Figures (ref), (ref), (ref) and (ref)), Figures (ref) shows that our proposed MoM-based estimator of $\alpha_{1}$ and $\alpha_{100}$ generally have similar $\sqrt{n} \times \mathrm{bias}$, variance, and mean squared error to the estimator proposed in bellec2025observable. Moreover, Figures (ref) and (ref) showcase the histograms and normal quantile-quantile plots of the $U$-statistic-based moment estimators and estimators of $\lambda_{\alpha}$, $\gamma_{\alpha}^{2}$, $\alpha_{1}$ and $\alpha_{100}$ by solving the moment system (ref). It is clear from these figures that the sampling distributions of both the $U$-statistic-based moment estimators and our estimators $\alpha_{1}$ and $\gamma_{\alpha}^{2}$ are close to the Gaussian distribution, further confirming by our theoretical results on the GAN property of our proposed MoM-based estimators. • When the design is Rademacher (Settings 3 & 4, Figures (ref), (ref), (ref) and (ref)), we observe similar results to those in the Gaussian settings. Therefore, under certain conditions on the design and the regression coefficients, our identification and estimation strategies based on Gaussian designs continue to be applicable -- demonstrating the universality of our proposed procedure. Interestingly, the debiased estimator also exhibits universality, which deserves a further theoretical investigation. It is worth noting that although Setting 4 concerns sparse regression coefficients, the number of non-zero coefficients is large and the values of non-zero coefficients are the same, so numerically our MoM-based estimators work still fine. In Appendix (ref), we showcase a different simulation setting, in which only one coordinate of the coefficients is non-zero; there the results show that the Gaussian-based identification and estimation strategies no longer produce consistent estimators. \end{itemize} \begin{figure}[htbp] \caption{Setting 1 of Section (ref) (Gaussian design and dense regression coefficients).} \end{figure} \begin{figure}[htbp] \caption{Setting 1 of Section (ref) (Gaussian design and dense regression coefficients): Sampling distributions of the moment estimators and the parameter estimators, over 500 Monte Carlos are displayed for the case of $n = 5000$.} \end{figure} \begin{figure}[htbp] \caption{Setting 3 of Section (ref) (Rademacher design and dense regression coefficients).} \end{figure} \begin{figure}[htbp] \caption{Setting 3 of Section (ref) (Rademacher design and dense regression coefficients): Sampling distributions of the moment estimators and the parameter estimators, over 500 Monte Carlos are displayed for the case of $n = 5000$.} \end{figure} \subsection{Mean estimation with missing data under MAR} For the problem of estimating the mean of a response $Y$ under Model (ref), we also consider several different settings. We first describe the common part of the DGP: the data is generated according to Model (ref), with $\mathbf{X}_{i} \overset{\rm i.i.d.}{\sim} \mathrm{N}_{p} (\bm{0}, \mathbf{I}_{p})$ for $i = 1, \cdots, n$. We take ${\mathbb{E}} [A | \mathbf{X}] = 0.1 + 0.9 \cdot \mathsf{expit} (\mathbf{X}^{\top} \bm{\alpha})$, as in celentano2023challenges. Here the true value of the target parameter $\psi$ is 0 because $\mathbf{X}$ has mean zero. We take the outcome noise $\varepsilon \sim \mathrm{N} (0, \sigma^{2})$ with $\sigma = 0.2$. The sample size varies from 100 to 10000. We compare our proposed MoM-based estimator $\widehat{\psi}$ with the two estimators proposed in celentano2023challenges, respectively termed as “Oracle ASCW” and “Empirical SCA”. For these alternative estimators, initial estimates based on Ridge are also used, With the tuning parameter $\lambda$ taking 50 values ranging from $[0.01, 10]$ with equal steps on a logarithmic scale. We refer readers to celentano2023challenges for the details of these two alternative estimators; here we simply implement the R code provided in the \href{https://github.com/mcelentano/Debiasing_for_missing_data}{GitHub repository} provided by the authors of celentano2023challenges. We also consider two different configurations for the regression coefficients $\bm{\alpha}$ and $\bm{\beta}$. In Setting 1 (Figures (ref) -- (ref)), we consider dense regression coefficients, where both $\bm{\alpha}$ and $\bm{\beta}$ are drawn i.i.d. coordinate-wise from $\mathrm{Uniform} ([- \sqrt{3 / p}, \sqrt{3 / p}])$. In Appendix (ref), we report results when changing the covariates distribution from Gaussian to $\mathbf{X}_{i, j} \overset{\rm i.i.d.}{\sim} \mathrm{Rad} (1 / 2)$ for $i = 1, \cdots, n$ and $j = 1, \cdots, p$. In Setting 2 (Figures (ref) -- (ref)), we consider sparse regression coefficients, where only $\alpha_{1}$ and $\beta_{1}$ are non-zero and are both equal to 1. Again, to save space, we defer the figures for Setting 2 to Appendix (ref). As can be seen from Figures (ref) and (ref), the MoM-based estimators generally have lower $\sqrt{n} \times \mathrm{bias}$, variance, and mean squared error than the two estimators proposed in celentano2023challenges, over a range of values of the tuning parameter $\lambda$, regardless of the configurations of the regression coefficients. Figures (ref) and (ref) display the histograms and normal quantile-quantile plots of the $U$-statistic-based moment estimators and the estimators of the target parameter $\psi$, together with $\lambda_{\alpha}, \gamma_{\alpha}^{2}, \gamma_{\alpha, \beta}$ (the byproducts of the system (ref)). It is clear from these figures that the sampling distributions of both the $U$-statistic-based moment estimators and of the estimators of the target parameter $\psi$, together with $\lambda_{\alpha}, \gamma_{\alpha}^{2}, \gamma_{\alpha, \beta}$ (the byproducts of the system (ref)), are close to the Gaussian distribution. \begin{figure}[htbp] \caption{Simulation results for Setting 1 (dense regression coefficients) in Section (ref). The two methods proposed in celentano2023challenges are plotted separately in two columns of the figure, with color gradients from blue to red representing the increasing value of the tuning parameter $\lambda$. The MoM-based estimators are plotted with white circles and dashed black lines.} \end{figure} \begin{figure}[htbp] \caption{Simulation results for Setting 1 (dense regression coefficients) in Section (ref): Sampling distributions of the moment estimators and the parameter estimators, over 500 Monte Carlos are displayed for the case of $n = 5000$.} \end{figure} \section{Discussion} This paper proposes Method-of-Moments (MoM) estimators of functionals of the regression coefficients in high-dimensional GLMs under the proportional asymptotic regime, in the most part assuming Gaussian designs with known population covariance matrix $\mathbf{\Sigma}$, following the current trend of literature. We demonstrate promising theoretical and numerical results about the MoM estimators. However, we emphasize that a more delicate comparison between our proposed approach and those based on debiasing bellec2025observable, celentano2023challenges is warranted in a future work: for instance, which estimator has better efficiency or which estimator performs better when $\mathbf{\Sigma}$ is unknown or when sufficient conditions (e.g. Assumption (ref)) for universality are violated with non-Gaussian designs. We end our article with the following discussion. \subsection{Statistical Inference} A problem that is worth pursuing in future work relates to the potentially non-Gaussian limiting distribution. When $\Vert \bm{\beta} \Vert$ is small, the corresponding $U$-statistic estimator may have positive probability of being negative and the truncation used in (ref) could lead to non-Gaussian limiting distributions of $\widehat{\gamma}_{\beta}^{2}$. It will be an interesting problem to investigate whether the bootstrap approach also delivers asymptotically valid inference in such a setting. \subsection{Unknown Population Covariance Matrix} In this paper, we mostly assume that $\mathbf{\Sigma}$ is known, an assumption ubiquitously adopted in the current literature (see Section (ref)). Progress can be made when $\mathbf{\Sigma}$ is unknown, but under specific conditions; see e.g. Section (ref) and Appendix (ref). It will be of greater interest to examine the case if we weaken the assumptions on $\mathbf{\Sigma}$ or $\mathbf{X}$ imposed in Section (ref). \subsection{Single-Index Models (SIM) and Model Misspecification} Compared to GLMs, isotonic Single-Index Models (SIMs) offer a more flexible modeling option as they allow the link function to be unknown and nonparametric. To the best of our knowledge, bellec2025observable is among the first to systematically study inference of functionals of regression coefficients in SIMs under proportional asymptotics. It is interesting to extend our framework to single-index models but it requires estimating the link function and hence higher-order moments, a complication that we plan to address in a separate paper. Another direction is to investigate if similar results for ATE still hold under the assumptions of su2023estimated, where only the propensity score is modeled by high-dimensional logistic regression but the outcome regression model is arbitrary. Finally, whenever parametric models such as GLMs are used, model misspecification is of potential concerns. The semiparametric partially nonlinear regression parameters considered in vansteelandt2022assumption, which reduce to functionals of regression coefficients when GLMs are correctly specified, can be another promising future avenue for our work. \section*{Acknowledgments} The authors would like to thank Subhabrata Sen, Pragya Sur, and Nicolas Verzelen for helpful discussions. The computations in this paper were run on the Siyuan-1 and $\pi$-2.0 cluster supported by the Center for High Performance Computing at Shanghai Jiao Tong University. Xingyu Chen and Lin Liu were supported by NSFC Grant No. 12471274. Lin Liu was also partly supported by NSFC Grant No. 12090024. \putbib[Master.bib]