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

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

115,224 characters

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



\title{Method-of-Moments Inference for GLMs and Doubly Robust Functionals under Proportional Asymptotics}

\author[1]{Xingyu Chen \thanks{E-mail: \texttt{[email removed]}}}
\author[1, 2]{Lin Liu \thanks{E-mail: \texttt{[email removed]}}}
\author[3]{Rajarshi Mukherjee \thanks{E-mail: \texttt{[email removed]} R. Mukherjee is supported by NSF CAREER 8529216-01.}}

\affil[1]{School of Mathematical Sciences, CMA-Shanghai, Shanghai Jiao Tong University}
\affil[2]{Institute of Natural Sciences, MOE-LSC, SJTU-Yale Joint Center for Biostatistics and Data Science, Shanghai Jiao Tong University}
\affil[3]{Department of Biostatistics, Harvard T. H. Chan School of Public Health}

\date{\today}

\maketitle
\begin{abstract}
In 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.
\end{abstract}


\allowdisplaybreaks

\begin{bibunit}[plainnat]

\section{Introduction}
\label{sec:intro}

Statistical inference in Generalized Linear Models (GLMs) \citep{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 \citep{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 \citep{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 \citep{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 \citep{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 \citep{negahban2012unified}. However, even the mere existence of such a consistent estimator relies on \textit{a priori} unknown low-dimensional (such as sparsity) assumptions in respective GLMs \citep{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 \citep{bean2013optimal, el2013asymptotic, donoho2016high, lei2018asymptotics}\footnote{In the econometrics literature, similar problems have also been studied in (partially) linear models \citep{cattaneo2018inference, cattaneo2019two} under the name ``many-regressor asymptotics'' or in settings with many weak instrumental variables \citep{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 \citep{verzelen2012minimax, barbier2019optimal}, recent research has coined it as the ``inconsistency regime'' \citep{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 \citep{bellec2023debiasing} and subsequent demonstration of a universality principle \citep{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 \textit{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{as:link} 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*)]
\item We propose moments-based identification strategies (for the precise meaning of identification, we refer readers to Lemma~\ref{lem:glm mean zero} 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.

\item 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.

\item 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.

\item 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.

\item 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}
\label{sec:review}
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 \citep{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 \citep{jankova2018semiparametric} and requires the consistent estimation of ultra high-dimensional (when the dimension $p$ is \textit{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 \citep{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 \textit{proportional} to the sample size $n$) \citep{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 \citet{bayati2011lasso, stojnic2013framework, thrampoulidis2018precise, miolane2021distribution, celentano2023lasso} and references therein) and then beyond Gaussian \citep{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 \citep{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 \citep{bellec2025observable}\footnote{It is noteworthy that \citet{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{sec:discussion}.}. Another work related to ours is \citet{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 \citet{sawaya2023moment} for a precise statement) and (2) the covariates $\mathbf{X}$ having zero mean, \citet{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 \cite{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 \cite{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 \citep{celentano2023challenges} or for specific functionals with $p < n$ \citep{yadlowsky2022explaining, jiang2025new}. Since the ultra-high-dimensional regime under sparsity has been heavily studied \citep{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 \citep{robins2008higher, van2021higher, kennedy2024minimax, bonvini2024doubly, breunig2019simple, breunig2024adaptive}. Also see Remark~\ref{rem:hoif} of Section~\ref{sec:unknown} 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 \citep{bellec2023debiasing, bellec2025observable}, in particular when $p > n$. Indeed, \citet{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 \citep{verzelen2018adaptive}. Its parallels in the proportional asymptotic regime without sparsity assumptions are quite sporadic, and we are only aware of \citet{li2023spectrum}, and to some extent \citet{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{sec:unknown}, 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 \citet{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 \citep{robins2008higher, robins2016technical}; we will discuss their similarity and difference further in Remark~\ref{rem:hoif}. 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{sec:glm} (knowing $\mathbf{\Sigma}$) and \ref{sec:unknown} (not knowing $\mathbf{\Sigma}$), we present our results on inference in GLMs followed by its applications in observational studies collected in Section~\ref{sec:obs}. Subsequently, Section~\ref{sec:sims} validates the theoretical results via numerical experiments. Our article ends by discussing some open problems in Section~\ref{sec:discussion}. 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}
\label{sec:glm}

In this section, we illustrate our main idea under the following stylized GLM. Suppose that we observe
\begin{equation}
\label{GLM} \tag{$\mathsf{GLM}$}
\begin{split}
& (Y_i, \mathbf{X}_i)_{i = 1}^{n} \stackrel{\rm i.i.d.}{\sim}\mathbb{P}_{\bm{\beta}}, \text{ 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 \textit{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}
\label{key parameters}
\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. \citet{guo2019optimal, song2024hede} and references therein. Moreover, as we will see while studying problems of estimating functionals of interest in observational studies in Section~\ref{sec:obs} 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{app:clt}).

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{tab:glossary}, together with where the assumptions are imposed throughout this paper.

\begin{table}[htbp]
\centering
\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{sec:zero mean} and \ref{sec:unknown} \\
$\mathsf{U}_{0}$ & Universality Conditions and Knowing $\bm{\mu} = \bm{0}$ & Sections~\ref{sec:zero mean} and \ref{sec:unknown} \\
$\mathsf{G}$ & Gaussian Design with Unknown $\bm{\mu}$ & Except Sections~\ref{sec:zero mean} and \ref{sec:unknown} \\
$\mathsf{U}$ & Universality Conditions with Unknown $\bm{\mu}$ & Except Sections~\ref{sec:zero mean} and \ref{sec:unknown} \\
\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}$.}
\label{tab:glossary}
\end{table}

First, we state the following ``global'' assumptions imposed on $\bm{\mu}$,  $\mathbf{\Sigma}$, and the link functions $\phi$.

\begin{customas}{$\mathsf{D}$}
\label{as:Sigma}
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}$}
\label{as:link}
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}
\item[(1)] $\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$.
\item[(2)] $\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}
\label{rem:link}
Assumption~\ref{as:link} 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{as:link}~(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{app:mu}.
\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}$}
\label{as:proportion}
There exists $\delta \in (0, \infty)$ such that $\lim_{n \rightarrow \infty} p / n \rightarrow \delta$.
\end{customas}

\begin{remark}
\label{rem:assumptions}
To shorten the exposition, we always assume Assumptions~\ref{as:Sigma},~\ref{as:link}, and~\ref{as:proportion} 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{rem:boundary} later.
\begin{customas}{$\mathsf{B}$}\leavevmode
\label{as:bounded}
\begin{itemize}
\item[(1)] There exists a universal constants $0 < \bar{B} < \infty$ such that $\Vert \bm{\beta} \Vert \leq \bar{B}$;
\item[(2)] 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
\label{as:var-cov}
\begin{itemize}
\item[(1)] 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}$;
\item[(2)] 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{as:var-cov}(2) is mainly made to ease exposition when we establish the CAN property of our proposed estimator (see e.g. Proposition~\ref{prop:glm_clt} and Theorem~\ref{thm:GLM CLT}). 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{sec:obs}, 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}
\label{sec:zero mean}

To gather intuition for our method, it is instructive first to consider the following assumptions.

\begin{customas}{$\mathsf{G_0}$}
\label{as:normal mean zero}
$\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}
\label{lem:glm mean zero}
Under Model~\ref{GLM}, Assumptions~\ref{as:normal mean zero},~\ref{as:bounded}{\rm(1)} and~\ref{as:var-cov}{\rm(1)}, the following hold:
\begin{enumerate}
\item[\emph{(1)}] 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}
\label{original chain}
\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}
\label{GLM mean zero chain}
\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}, \label{mm} \\
& 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}, \label{mm2}
\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 \eqref{mm} as $\Psi_{\mathsf{GLM}_{0}, \beta}: \gamma_{\beta}^{2} \mapsto m_{\mathbf{X} Y, 2}$ and the map induced by \eqref{GLM mean zero chain} 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}})$.

\item[\emph{(2)}] 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 \eqref{GLM chain} uniquely determines the value of $(\gamma_{\beta}^{2}, \bm{\beta}^{\top} \bm{\upsilon})$.
\end{enumerate}
\end{lemma}

We now unpack Lemma~\ref{lem:glm mean zero}, with its proof deferred to the Appendix. Based on the Gaussian design and Stein's lemma, the moment equations \eqref{original chain} 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 \eqref{mm}, that maps the quadratic form $\gamma_{\beta}^{2}$ to the moment $m_{\mathbf{X} Y, 2}$. Assumption~\ref{as:link} on the link function and Assumption~\ref{as:bounded} 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 \eqref{mm2}, $\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{lem:glm mean zero}, together with all the other identification results under Gaussian designs in this paper, do not require Assumption~\ref{as:proportion}. 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. \cite{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]
\label{def:id}
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{GLM}) 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
\label{as:beta mean zero}
\begin{enumerate}[label = (\arabic*)]
\item $\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$;
\item $\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{as:beta mean zero}(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{lem:glm mean zero}, without assuming that $\mathbf{X}$ is Gaussian.

\begin{lemma}
\label{lem:glm mean zero, universality}
Under the same assumptions as in Lemma~\ref{lem:glm mean zero}, except with Assumption~\ref{as:normal mean zero} replaced by Assumption~\ref{as:beta mean zero}, the system of moment equations appeared in Lemma~\ref{lem:glm mean zero} 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{app:universality}. Taken together, the above results motivate the following MoM-based estimator of $\gamma_{\beta}^{2}$ and $\bm{\beta}^{\top} \bm{\upsilon}$:
\begin{equation}
\label{key estimator}
\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{app:universality}). We next summarize the above reasoning as the following proposition.

\begin{proposition}
\label{prop:glm_rate}
Under the Assumptions of Lemma~\ref{lem:glm mean zero} or Lemma~\ref{lem:glm mean zero, universality}, 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}
\label{prop:glm_clt}
Under Model~\ref{GLM}, Assumptions~\ref{as:normal mean zero},~\ref{as:bounded} and~\ref{as:var-cov}, 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}
\label{rem:boundary}
It is noteworthy that Assumption~\ref{as:bounded}(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}
\label{rem:clt}
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 \citet{bhattacharya1992class} (see Appendix~\ref{app:CLT} 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}
\label{rem:clt condition}
In Proposition~\ref{prop:glm_clt} (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{app:CLT}, 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{as:beta mean zero}, 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{as:var-cov}. 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{lem:clt} and Proposition~\ref{prop:general clt} in Appendix~\ref{app:CLT}).

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 \citet{liu2024assumption}), consistent variance estimators based on tweaking the nonparametric bootstrap have been developed and can be used to conduct inference; also see Section~\ref{sec:bootstrap} for a brief discussion and Appendix~\ref{app:bootstrap} for its finite-sample performance.
\end{remark}

\begin{remark}
\label{rem:debunk}
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 \citet{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}
\label{sec:mu}

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}$}
\label{as:normal}
$\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{lem:glm mean zero} to the case of Gaussian designs with unknown mean $\bm{\mu}$.
\begin{lemma}
\label{lem:mu}
Under Model~\ref{GLM}, Assumptions~\ref{as:normal},~\ref{as:bounded}{\rm(1)}, and~\ref{as:var-cov}{\rm(1)}:
\begin{enumerate}
\item[\emph{(1)}] The following system of moment equations holds:
\begin{subequations}
\label{GLM chain}
\begin{align}
& m_{Y} \coloneqq {\mathbb{E}} [Y] = \mathrm{f}_{0} (\textcolor{black}{\lambda_{\beta}}, \textcolor{black}{\gamma_{\beta}^{2}}), \label{GLM m1} \\
& 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}, \label{GLM m2} \\
& 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}}, \label{GLM m3} \\
& 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}}, \label{GLM m4} \\
& 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}, \label{GLM m5} \\
& 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}. \label{GLM m6}
\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*}

\item[\emph{(2)}] Further, the first four equations of \eqref{GLM chain}, denoted as $\Psi_{\mathsf{GLM}, [4]}$, can be reduced to
\begin{subequations}
\label{GLM chain reduced}
\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 \eqref{GLM chain reduced} 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}
\label{rem:difficult}
The proof of Lemma~\ref{lem:mu} can be found in Appendix~\ref{app:mu} and Appendix~\ref{app:non-zero}. As mentioned in the beginning of this section, the system of moment equations in \eqref{GLM chain reduced} is from ${\mathbb{R}}^{2}$ to ${\mathbb{R}}^{2}$. As a result, showing that \eqref{GLM} is a global diffeomorphism is nontrivial, which involves the use of Hadamard global inverse function theorem. For details, see Lemma~\ref{lem:inverse} and Appendix~\ref{app:non-zero}.
\end{remark}

To establish universality of Lemma~\ref{lem:mu} beyond Gaussian designs, we need to first generalize Assumption~\ref{as:beta mean zero} to Assumption~\ref{as:beta} below.

\begin{customas}{$\mathsf{U}$}\leavevmode
\label{as:beta}
\begin{enumerate}[label = (\arabic*)]
\item $\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$;
\item $\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{lem:glm mean zero, universality} can then be generalized as follows, the proof of which can be found in Appendix~\ref{app:universality}.

\begin{lemma}
\label{lem:mu, universality}
Under the same assumptions as in Lemma~\ref{lem:mu}, except with Assumption~\ref{as:normal} replaced by Assumption~\ref{as:beta}, the systems of moment equations appeared in Lemma~\ref{lem:mu} 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{as:beta}\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{as:beta}\rm{(2)} is violated -- as demonstrated via the numerical experiments related to Figures~\ref{fig:GLM, sparse-converge, rademacher, univ} and~\ref{fig:GLM, hist, sparse-converge, rademacher, univ}; see Section~\ref{sec:sims} for more details.

\begin{remark}
\label{rem:zero mean}
Compared to the special case of knowing $\bm{\mu} \equiv \bm{0}$ in Section~\ref{sec:zero mean}, it requires extra moment equations \eqref{GLM m1} to \eqref{GLM m3}  to identify $\gamma_{\beta}^{2}$, together with the linear form $\lambda_{\beta}$. When it is known that $\bm{\mu} = \bm{0}$, \eqref{GLM m1}, \eqref{GLM m2}, and \eqref{GLM m3}, respectively, reduce to constants $1 / 2$, $0$, and $0$.
\end{remark}

\begin{remark}
\label{rem:non-gaussian}
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}
\label{joint estimator}
\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}
\label{GLM moments estimators}
\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 \eqref{GLM m5} to \eqref{GLM m6}:
\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}
\label{thm:GLM}
Under the Assumptions of Lemma~\ref{lem:mu} or Lemma~\ref{lem:mu, universality}, 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{prop:glm_clt} for the CAN property of our proposed estimators to the case where $\bm{\mu}$ is unknown and possibly non-zero.


\begin{theorem}
\label{thm:GLM CLT}
Under Model~\ref{GLM}, Assumptions~\ref{as:normal},~\ref{as:bounded} and~\ref{as:var-cov}, 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{app:GLM functional CLT}, 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{thm:GLM CLT} by the bootstrap method developed in \citet{liu2024assumption}. The performance of bootstrap variance estimators will be assessed in Appendix~\ref{app:bootstrap}.

\begin{remark}[On the scalings between $p$ and $n$]
\label{rem:p}
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 \citet{bhattacharya1992class}.
\end{remark}

\begin{remark}[Data-driven Approximate Message Passing Schemes for GLMs]
\label{rem:amp}
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{sec:obs}. For logistic regression, i.e. Model \eqref{GLM} with $\phi \equiv \mathsf{expit}$, \citet{sur2019modernb} (and in a more general form, \citet{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}
\label{fixed point 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 \citet{sur2019modernb} leveraging and \texttt{SLOE} by \citet{yadlowsky2021sloe} using a reparameterization of \eqref{fixed point equation} 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 \citet{sawaya2023moment} try to solve. However, as mentioned in the Introduction, the approach in \citet{sawaya2023moment} only applies to asymmetric link functions because it does not rely on $\mathbf{\Sigma}$.
\end{remark}



\section{The Case of Unknown Population Covariance}
\label{sec:unknown}

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{app:unknown complete} 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}}
\label{sec:unknown well-posed}

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}$ \citep{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 \eqref{key estimator}, except that $\widehat{m}_{\mathbf{X} Y, 2}$ is replaced by
\begin{equation}
\label{alternative}
\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{app:unknown GLMs}}, $\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{prop:glm_rate}, the proof of which can be found in Appendix~\ref{app:unknown GLMs}.

\begin{proposition}
\label{prop:unknown}
Under Assumptions of Lemma~\ref{lem:glm mean zero}, 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{lem:Wishart} for more details. It is worth noting that Proposition~\ref{prop:unknown} can also be generalized beyond Gaussian designs by using universality-type results on (inverse) sample Gram matrices, e.g. related advances in \citet{derezinski2021sparse} and \citet{derezinski2022unbiased}.

\begin{remark}
\label{rem:unknown linear}
Under linear models, one can in fact adapt the estimators of \citet{guo2022moderate} for unknown $\mathbf{\Sigma}$ by leveraging the linearity structure; for details, see Appendix~\ref{app:unknown linear}.

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}}
\label{sec:unknown ill-posed}

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 \citet{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 \citep{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{app:ill-posed}. 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{app:ill-posed}. 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}
\label{prop:ill posed}
Under Assumptions of Lemma~\ref{lem:glm mean zero} or Lemma~\ref{lem:glm mean zero, universality}, 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}
\label{rem:hoif}
$\widehat{m}_{\mathbf{X} Y, 2}$ defined in \eqref{alternative} 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 \citet{liu2017semiparametric}. In \citet{liu2017semiparametric}, however, they did not make distributional assumptions on $\mathbf{X}$ such as Assumption~\ref{as:normal mean zero} or even Assumption~\ref{as:normal}, 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 \citep{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 \citet{kong2018estimating} to avoid estimating $\mathbf{\Sigma}^{-1}$. \citet{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{sec:mu} and Appendix~\ref{app:unknown complete}), making the computation a much more challenging problem that we plan to address in a separate paper. We mention in passing that \citet{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}
\label{sec:obs}

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{sec:unknown}. The universality of the proposed procedure can be established by analogous arguments in the proof of Lemma~\ref{lem:mu, universality}, 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 \citep{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 \citet{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 \citep{liu2024assumption}) that have gained popularity in recent literature on conditional independence testing \citep{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{sec:MAR} with a minor modification. The approaches considered in Section~\ref{sec:unknown} 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}
\label{sec:CE}

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}
\label{CE} \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{CE}, and the derivation can be found in Appendix~\ref{app:identification}.

\begin{lemma}
\label{lem:CE}
Under Model~\ref{CE}, Assumptions~\ref{as:normal},~\ref{as:bounded}{\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and~\ref{as:var-cov}{\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}, A$, the following system of moment equations holds:
\begin{subequations}
\label{CE chain}
\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{lem:mu}.

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{app:CE}. As a direct corollary of Lemma~\ref{lem:CE}, 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 \eqref{GLM moments estimators}.

\begin{theorem}
\label{thm:CE}
Under Model~\ref{CE}, Assumptions~\ref{as:normal},~\ref{as:bounded}{\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and~\ref{as:var-cov}{\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{CE} with $\varepsilon$ being independent from $\mathbf{X}, A$, Assumptions~\ref{as:normal},~\ref{as:bounded} on both $\bm{\alpha}$ and $\bm{\beta}$, and~\ref{as:var-cov} 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, \cite{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{sec:obs}, including Theorem~\ref{thm:CE}, can be found in Appendix~\ref{app:CLT}. 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{thm:GLM CLT}.



\subsection{Mean Estimation with Missing Data under Missing-At-Random (MAR)}
\label{sec: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}
\label{MAR} \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 \citet{celentano2023challenges}, we assume that only $\mathbf{\Sigma}$ is known. Under Model~\ref{MAR},
\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{MAR}, and the derivation can again be found in Appendix~\ref{app:identification}.
\begin{lemma}
\label{lem:MAR}
Under Model~\ref{MAR}, Assumptions~\ref{as:normal},~\ref{as:bounded}{\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and~\ref{as:var-cov}{\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, the following system of moment equations holds:
\begin{subequations}
\label{MAR chain}
\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}}, \label{mar e} \\
& 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}}, \label{mar f}
\end{align}
\end{subequations}
where the definition of $\mathrm{f}_{k}$ appears in the statement of Lemma~\ref{lem:mu}.

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}
\label{alt}\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 \eqref{alt} with either \eqref{mar e} or \eqref{mar f} 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 \eqref{MAR chain} 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{app:MAR}. As a direct corollary of Lemma~\ref{lem:MAR}, 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 \eqref{GLM moments estimators}.

\begin{theorem}
\label{thm:MAR}
Under Model~\ref{MAR}, Assumptions~\ref{as:normal},~\ref{as:bounded}{\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and~\ref{as:var-cov}{\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{MAR}, Assumptions~\ref{as:normal},~\ref{as:bounded} on both $\bm{\alpha}$ and $\bm{\beta}$, and~\ref{as:var-cov} 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{MAR}, as we mentioned in the Introduction, \citet{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 \citet{celentano2023challenges} later in Section~\ref{sec:sims mar}.



\subsection{Estimation of the Generalized Covariance Measure (GCM)}
\label{sec:GCM}

In this section, we assume that ${\mathbb{P}}$ is defined and parameterized as follows:
\begin{equation}
\label{GCM} \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)$. \citet{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{GCM}, and the derivation can again be found in Appendix~\ref{app:identification}.

\begin{lemma}
\label{lem:GCM}
Under Model~\ref{GCM}, Assumptions~\ref{as:normal},~\ref{as:bounded}{\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$, and~\ref{as:var-cov}{\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, the following system of moment equations holds:
\begin{subequations}
\label{GCM chain}
\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} \label{mm3} \\
& + 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{lem:mu}.

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{app:GCM}. As a direct corollary of Lemma~\ref{lem:GCM}, 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 \eqref{GLM moments estimators}.

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
\label{as:cov}
\begin{itemize}
\item[(1)] We assume that $\Vert \sigma_{A, Y} (\cdot) \Vert_{2}$ is bounded;
\item[(2)] 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}
\label{thm:GCM}
Under Model~\ref{GCM}, Assumptions~\ref{as:normal},~\ref{as:bounded}{\rm(1)} on both $\bm{\alpha}$ and $\bm{\beta}$,~\ref{as:var-cov}{\rm(1)} on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, and~\ref{as:cov}\rm{(1)}, we have
\begin{align*}
\sqrt{n} (\widehat{\psi} - \psi) = O_{{\mathbb{P}}} (1).
\end{align*}
Alternatively, under Model~\ref{GCM}, Assumptions~\ref{as:normal},~\ref{as:bounded} on both $\bm{\alpha}$ and $\bm{\beta}$,~\ref{as:var-cov} on both $A | \mathbf{X}$ and $Y | \mathbf{X}$, and~\ref{as:cov}, 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{as:cov}(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{sec:CE}, 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{sec:MAR}, we assume that $Y \protect\mathpalette{\protect\independenT}{\perp} A | \mathbf{X}$.



\section{Numerical Experiments}
\label{sec:sims}

In this section, we verify our theory via extensive numerical experiments. We consider two types of problems. The first problem, studied in Section~\ref{sec:sims glms}, 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{sec:sims mar}, is on estimating the mean of a response $Y$ subject to missingness under Model~\ref{MAR}. 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{sec:sims glms} and 1.25 in Section~\ref{sec:sims mar}. Appendix~\ref{app:sim} 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{app:sims unknown} we conduct simulations to evaluate the performance of our proposed estimators in Section~\ref{sec:unknown well-posed} when knowing neither $\bm{\mu}$ nor $\mathbf{\Sigma}$. The bootstrap variance estimators mentioned right after Theorem~\ref{thm:GLM CLT} are assessed in Appendix~\ref{app:bootstrap}.

For all the histograms and normal quantile-quantile plots reported in this section and in Appendix~\ref{app:sim}, 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}
\label{sec:sims 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}
\item $\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.
\item $\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{app:sim fig}.
\item 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 \citet{bellec2025observable}.
\end{itemize}

In each setting, for the problem of estimating certain single coordinate, we compare our MoM-based estimators against the methods of \citet{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 \citet{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}
\item When the design is Gaussian (Settings 1 \& 2, Figures~\ref{fig:GLM, dense-converge, gaussian，1.2},~\ref{fig:GLM, dense-converge, gaussian，1.2, hist},~\ref{fig:GLM, sparse-converge, gaussian，1.2} and~\ref{fig:GLM, sparse-converge, gaussian，1.2, hist}), Figures~\ref{fig:GLM, dense-converge, gaussian，1.2} 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 \citet{bellec2025observable}. Moreover, Figures~\ref{fig:GLM, dense-converge, gaussian，1.2, hist} and~\ref{fig:GLM, sparse-converge, gaussian，1.2, hist} 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 \eqref{GLM chain}. 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.

\item When the design is Rademacher (Settings 3 \& 4, Figures~\ref{fig:GLM, dense-converge, rademacher，1.2},~\ref{fig:GLM, dense-converge, rademacher，1.2, hist},~\ref{fig:GLM, sparse-converge, rademacher，1.2} and~\ref{fig:GLM, sparse-converge, rademacher，1.2, hist}), 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{app:sim univ}, 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]
\centering
\includegraphics[width = 0.65\textwidth]{./Simulations_version_1/simu500-rootn_bias_glm_mom_sparse_0_one_0_Ra_0_p_over_n_1.2.pdf}
\caption{Setting 1 of Section~\ref{sec:sims glms} (Gaussian design and dense regression coefficients).}
\label{fig:GLM, dense-converge, gaussian，1.2}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width = 0.65\textwidth]{./Simulations_version_1/glm_mom_norm_sparse_0_one_0_Ra_0_p_over_n_1.2.pdf}
\caption{Setting 1 of Section~\ref{sec:sims glms} (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$.}
\label{fig:GLM, dense-converge, gaussian，1.2, hist}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width = 0.65\textwidth]{./Simulations_version_1/simu500-rootn_bias_glm_mom_sparse_0_one_0_Ra_1_p_over_n_1.2.pdf}
\caption{Setting 3 of Section~\ref{sec:sims glms} (Rademacher design and dense regression coefficients).}
\label{fig:GLM, dense-converge, rademacher，1.2}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width = 0.65\textwidth]{./Simulations_version_1/glm_mom_norm_sparse_0_one_0_Ra_1_p_over_n_1.2.pdf}
\caption{Setting 3 of Section~\ref{sec:sims glms} (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$.}
\label{fig:GLM, dense-converge, rademacher，1.2, hist}
\end{figure}

\subsection{Mean estimation with missing data under MAR}
\label{sec:sims mar}

For the problem of estimating the mean of a response $Y$ under Model~\ref{MAR}, we also consider several different settings. We first describe the common part of the DGP: the data is generated according to Model~\ref{MAR}, 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 \citet{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 \citet{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 \citet{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 \citet{celentano2023challenges}.

We also consider two different configurations for the regression coefficients $\bm{\alpha}$ and $\bm{\beta}$. In Setting 1 (Figures~\ref{fig:MAR, dense, gaussian} --~\ref{fig:MAR, parameters, dense, gaussian}), 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{app:sim mar}, 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{fig:MAR, sparse, gaussian} --~\ref{fig:MAR, parameters, sparse, gaussian}), 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{app:sim mar}.

As can be seen from Figures~\ref{fig:MAR, dense, gaussian} and~\ref{fig:MAR, sparse, gaussian}, the MoM-based estimators generally have lower $\sqrt{n} \times \mathrm{bias}$, variance, and mean squared error than the two estimators proposed in \citet{celentano2023challenges}, over a range of values of the tuning parameter $\lambda$, regardless of the configurations of the regression coefficients. Figures~\ref{fig:MAR, parameters, dense, gaussian} and~\ref{fig:MAR, parameters, sparse, gaussian} 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 \eqref{MAR chain}). 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 \eqref{MAR chain}), are close to the Gaussian distribution.

\begin{figure}[htbp]
\centering
\includegraphics[width = 0.65\textwidth]{Simulations_version_1/glm_mar_lambda_sparse_0_Ra_0_p_over_n_1.25.pdf}
\caption{Simulation results for Setting 1 (dense regression coefficients) in Section~\ref{sec:sims mar}. The two methods proposed in \citet{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.}
\label{fig:MAR, dense, gaussian}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width = 0.65\textwidth]{Simulations_version_1/glm_mar_para_norm_sparse_0_Ra_0_p_over_n_1.25.pdf}
\caption{Simulation results for Setting 1 (dense regression coefficients) in Section~\ref{sec:sims mar}: Sampling distributions of the moment estimators and the parameter estimators, over 500 Monte Carlos are displayed for the case of $n = 5000$.}
\label{fig:MAR, parameters, dense, gaussian}
\end{figure}

\section{Discussion}
\label{sec: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 \citep{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{as:beta}) for universality are violated with non-Gaussian designs. We end our article with the following discussion.

\subsection{Statistical Inference}
\label{sec:bootstrap}


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 \eqref{key estimator} 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{sec:review}). Progress can be made when $\mathbf{\Sigma}$ is unknown, but under specific conditions; see e.g. Section~\ref{sec:unknown} and Appendix~\ref{app:unknown}. 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{sec:unknown}.

\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, \citet{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 \citet{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 \citet{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]
\end{bibunit}

\newpage