High-Dimensional Inference for Generalized Linear Models
with Hidden Confounding

Statistical inferences for high-dimensional regression models have received a growing interest due to the increasing number of complex data sets with high-dimensional covariates that are collected across many different scientific disciplines ranging from genomics to social science to econometrics (Peng et al., 2010; Belloni et al., 2012; Fan et al., 2014; Bühlmann et al., 2014). One of the most popularly used high-dimensional regression methods is the lasso linear regression, which assumes that the underlying regression coefficients are sparse (Tibshirani, 1996). However, the lasso penalty introduces non-negligible bias that renders high-dimensional statistical inference challenging (Wainwright, 2019). To address this challenge, Zhang and Zhang (2014) and van de Geer et al. (2014) proposed the debiasing method to construct confidence intervals for the lasso regression coefficients: the main idea is to first obtain the lasso estimator, and then correct for the bias of the lasso estimator using a low-dimensional projection method. We refer the reader to Javanmard and Montanari (2014), Belloni et al. (2014), Ning and Liu (2017), Chernozhukov et al. (2018), among others, for detailed discussions of the debiasing approach, and also Wainwright (2019) for an overview of high-dimensional statistical inference.

The aforementioned studies were established under the assumption that there are no unmeasured confounders that are associated with both the response and covariates. However, this assumption is often violated in observational studies. For instance, in genetic studies, the effect of certain segments of DNA on the gene expression may be confounded by population structure and microarray expression artifacts (Listgarten et al., 2010). Another example is in healthcare studies where the effect of nutrients intake on the risk of cancer may be confounded by physical wellness, social class, and behavioral factors (Fewell et al., 2007). Without adjusting for the unmeasured confounders, the resulting inferences from the standard debiasing methods could be biased and consequently, lead to spurious scientific discoveries.

Various methods have been proposed to perform valid statistical inferences for regression parameters in the presence of hidden confounders. One commonly used approach is instrumental variables regression, which typically requires domain knowledge to identify valid instrumental variables, making it challenging for high-dimensional applications (Kang et al., 2016; Burgess et al., 2017; Guo et al., 2018; Windmeijer et al., 2019). Recently, under a relaxed linear structural equation modeling framework, Guo et al. (2022) proposed a deconfounding approach for statistical inference for individual regression coefficients. Fan et al. (2023) developed a factor-adjusted debiasing method and conducted statistical estimation and inference for the coefficient vector. However, these two methods rely on the linearity assumption between the response and covariates, which may not hold when the response is not continuous but categorical (binary, count, etc.). Such categorical data commonly occurs in genetic, biomedical, and social science applications, and thus alternative methods are needed to perform valid statistical inferences in these cases.

To bridge the gap in the existing literature, we propose a novel framework to perform statistical inference in the context of high-dimensional generalized linear models with hidden confounding. The main idea of our proposed method involves estimating the unmeasured confounders using high-dimensional factor analysis techniques and then applying the debiasing method on the regression coefficient of interest, with the estimated unmeasured confounders treated as surrogate variables. This method does not rely on the linear model assumption or any specific model form and it is applicable to more general models beyond the generalized linear models, such as graphical models (Ren et al., 2015; Zhu et al., 2020) and additive hazards models (Lin and Ying, 1994; Lin and Lv, 2013).

Theoretically, we show that under some mild scaling conditions, the estimation errors of proposed estimators achieve comparable rates as that of the $\ell_{1}$-penalized generalized linear model without unmeasured confounders. We further show that the debiased estimator for the coefficient of interest is asymptotically Gaussian after adjusting for the unmeasured confounders, which results in a valid statistical inference for high-dimensional generalized linear models with hidden confounding. It is worth highlighting that when using a factor model to relate covariates and unmeasured confounders, we make more general assumptions on the random noise compared to existing works. Specifically, we allow the random noise to be non-identically distributed. This represents a significant improvement over the majority of previous works, which assume that the random noise follows an independent and identical distribution (Guo et al., 2022; Fan et al., 2023). Furthermore, unlike existing methods under the linear framework (Guo et al., 2022; Fan et al., 2023), generalized linear models pose new challenges as many transformation and decomposition techniques commonly used for the linear models are not directly applicable. Consequently, more general and refined intermediate results are needed as “building blocks” in the proofs (see Remark 8 in Appendix D.1 for details).

Our paper is organized and structured as follows. In Section 2, we introduce the model setup and provide a comprehensive discussion of related literature. Section 3 presents our two-step approach for estimating the parameter of interest while adjusting for the effect of unmeasured confounders. Section 4 establishes the theoretical properties of the parameter including the estimation consistency as well as the asymptotic normality for the estimator. In Section 5, we demonstrate the performance of our proposed method and the validity of theoretical results via extensive simulation studies. We also provide an application to genetic data containing gene expression quantifications and stimulation statuses in mouse bone marrow-derived dendritic cells, where we identify significant gene expressions under stimulations, which are consistent with the experimental findings in genetic studies (Nemetz et al., 1999; Ather and Poynter, 2018; Jang et al., 2018; Toyoshima et al., 2019). Lastly, in Section 6, we provide concluding remarks and outline potential future directions for this research.

In this section, we first set up a generalized linear model with hidden confounding and introduce a scientific application of our model framework. Then we will discuss related high-dimensional models in the existing literature.

In this section, we propose a novel framework to perform statistical inference for a parameter of interest in the context of high-dimensional generalized linear models with hidden confounding. In the proposed framework, we first estimate the unmeasured confounders using a factor analysis approach. Subsequently, the estimated unmeasured confounders are treated as surrogate variables for fitting a high-dimensional generalized linear model, and a debiased estimator is constructed to perform statistical inference.

Throughout this section, we assume that the observed data $\{y_{i},\bm{X}_{i}\}_{i=1,\dots,n}$ and the unmeasured confounders $\{\bm{U}_{i}\}_{i=1,\dots,n}$ are realizations of (1) and (2). Moreover, the random noise $\bm{E}_{i}=(E_{i1},\ldots,E_{ip})^{{\mathrm{\scriptscriptstyle T}}}$ has mean zero and variance $\bm{\Omega}_{i}=\mathbb{E}(\bm{E}_{i}\bm{E}_{i}^{{\mathrm{\scriptscriptstyle T}}}).$ Let $\bm{\Sigma}_{e}={\rm diag}(n^{-1}\sum_{i=1}^{n}\bm{\Omega}_{i})$, where ${\rm diag}(\bm{A})$ denotes a diagonal matrix by setting off-diagonal entries in $\bm{A}$ to zero. In the $p\times p$ diagonal matrix $\bm{\Sigma}_{e}$, we denote the $j$-th diagonal element to be $\sigma_{j}^{2}=n^{-1}\sum_{i=1}^{n}\tau_{i,jj}$, where $\tau_{i,jj}$ is the $(j,j)$ element of $\bm{\Omega}_{i}$. The model assumption on the random noise is general as it does not assume that the random noise $\bm{E}_{i}$ is identical nor does it require the covariance matrix $\bm{\Omega}_{i}$ to be diagonal. The detailed theoretical assumptions regarding the random noise will be presented in Section 4.

Estimation of the Unmeasured Confounders: In this work, we consider the dimension $K$ as pre-specified. In practice, there are various methods to estimate the dimension of unmeasured confounders such as scree plot (Cattell, 1966), cross-validation method (Owen and Wang, 2016), information criteria method (Bai and Ng, 2002), the eigenvalue ratio method (Lam and Yao, 2012; Ahn and Horenstein, 2013), among others. In the implementation of our proposed method, we recommend the parallel analysis approach because of its good finite sample performance, easy implementation, and popularity in scientific applications (Hayton et al., 2004; Costello and Osborne, 2005; Brown, 2015) – it has shown superior performances compared to many other methods in various empirical studies (Zwick and Velicer, 1986; Peres-Neto et al., 2005). Detailed discussions on the estimation of the dimension of unmeasured confounders are presented in Appendix C.

We employ the maximum likelihood estimation to estimate the unmeasured confounders under (2). Without loss of generality, let $\bar{\bm{U}}=n^{-1}\sum_{i=1}^{n}\bm{U}_{i}={0}$ and let $\bm{S}_{u}=n^{-1}\sum_{i=1}^{n}\bm{U}_{i}\bm{U}_{i}^{\mathrm{\scriptscriptstyle T}}$ be the sample variance of $\bm{U}$. Similarly, let $\bm{S}_{x}=n^{-1}\sum_{i=1}^{n}(\bm{X}_{i}-\bar{\bm{X}})(\bm{X}_{i}-\bar{\bm{X}})^{{\mathrm{\scriptscriptstyle T}}}$ be the sample variance of $\bm{X}$, where $\bar{\bm{X}}=n^{-1}\sum_{i=1}^{n}\bm{X}_{i}$. Given unmeasured confounders $\bm{U}_{1},\dots,\bm{U}_{n}$, define an approximation of population variance of $\bm{X}$ to be $\bm{\Sigma}_{x}=\bm{W}^{\mathrm{\scriptscriptstyle T}}\bm{S}_{u}\bm{W}+\bm{\Sigma}_{e}$. Note that $\bm{\Sigma}_{x}$ is not exactly the covariance matrix of $\bm{X}$ because we do not restrict $\bm{\Omega}_{i}$ to be diagonal and define $\bm{\Sigma}_{e}$ to be diagonal by setting the off-diagonal of $n^{-1}\sum_{i=1}^{n}\bm{\Omega}_{i}$ to be zero. Based on the factor model in (2), the maximum likelihood estimators of $\bm{W}$, $\bm{S}_{u}$ and ${\bm{\Sigma}}_{e}$ are obtained as follows:

$$(\hat{\bm{W}},\hat{\bm{S}}_{u},\hat{\bm{\Sigma}}_{e})=\underset{{\bm{W}},{\bm{S}}_{u},{\bm{\Sigma}}_{e}}{\operatorname{argmax}}\left\{-(2p)^{-1}\log|\bm{\Sigma}_{x}|-(2p)^{-1}\operatorname{tr}(\bm{S}_{x}\bm{\Sigma}_{x}^{-1})\right\}.$$

Computationally, we employ the Expectation-Maximization (EM) algorithm to obtain the maximum likelihood estimators as suggested in Bai and Li (2012) and Bai and Li (2016), where the authors proved the EM solutions are the stationary points for the likelihood function. Specifically, in this iterated EM algorithm, we use the principal components estimator as the initial estimator. Because principal component estimators are shown to be consistent estimators under similar model assumptions as ours (Fan et al., 2013; Wang and Fan, 2017), using the principal component estimators in initialization instead of using random initialization helps to improve algorithm efficiency and find more refined estimation results. At the $t$-th iteration, denote the estimators at this step to be $\bm{W}^{(t)}$ and $\bm{\Sigma}_{e}^{(t)}$. The EM algorithm updates the estimators to be

	$\displaystyle\{\bm{W}^{(t+1)}\}^{{\mathrm{\scriptscriptstyle T}}}$	$\displaystyle=$	$\displaystyle\left[n^{-1}\sum_{i=1}^{n}\mathbb{E}(\bm{X}_{i}\bm{U}_{i}^{{\mathrm{\scriptscriptstyle T}}}\mid\bm{X},\bm{W}^{(t)},\bm{\Sigma}_{e}^{(t)})\right]\left[n^{-1}\sum_{i=1}^{n}\mathbb{E}(\bm{U}_{i}\bm{U}_{i}^{{\mathrm{\scriptscriptstyle T}}}\mid\bm{X},\bm{W}^{(t)},\bm{\Sigma}_{e}^{(t)})\right]^{-1}$
	$\displaystyle\bm{\Sigma}_{e}^{(t+1)}$	$\displaystyle=$	$\displaystyle{\rm diag}(\bm{S}_{x}-\bm{W}^{(t+1)}\{\bm{W}^{(t)}\}^{{\mathrm{\scriptscriptstyle T}}}\{\bm{\Sigma}_{x}^{(t)}\}^{-1}\bm{S}_{x})$

where $\bm{\Sigma}_{x}^{t}=\{\bm{W}^{(t)}\}^{{\mathrm{\scriptscriptstyle T}}}\bm{W}^{(t)}+\bm{\Sigma}_{e}^{(t)}$,

	$\displaystyle n^{-1}\sum_{i=1}^{n}\mathbb{E}(\bm{X}_{i}\bm{U}_{i}^{{\mathrm{\scriptscriptstyle T}}}\mid\bm{X},\bm{W}^{(t)},\bm{\Sigma}_{e}^{(t)})$	$\displaystyle=$	$\displaystyle\bm{S}_{x}\{\bm{\Sigma}_{e}^{(t)}\}^{-1}\{\bm{W}^{(t)}\}^{{\mathrm{\scriptscriptstyle T}}},$
	$\displaystyle n^{-1}\sum_{i=1}^{n}\mathbb{E}(\bm{U}_{i}\bm{U}_{i}^{{\mathrm{\scriptscriptstyle T}}}\mid\bm{X},\bm{W}^{(t)},\bm{\Sigma}_{e}^{(t)})$	$\displaystyle=$	$\displaystyle\bm{W}^{(t)}\{\bm{\Sigma}_{e}^{(t)}\}^{-1}\bm{S}_{x}\{\bm{\Sigma}_{e}^{(t)}\}^{-1}\{\bm{W}^{(t)}\}^{{\mathrm{\scriptscriptstyle T}}}$
			$\displaystyle+\mathbb{I}_{K}-\bm{W}^{(t)}\{\bm{\Sigma}_{e}^{(t)}\}^{-1}\{\bm{W}^{(t)}\}^{{\mathrm{\scriptscriptstyle T}}}.$

The iterative steps stop when $\|\bm{W}^{(t+1)}-\bm{W}^{(t)}\|_{F}$ and $\|\bm{\Sigma}_{e}^{(t+1)}-\bm{\Sigma}_{e}^{(t)}\|_{F}$ are less than certain tolerance value. Let $\bm{W}^{{\dagger}}$ and $\bm{\Sigma}_{e}^{{\dagger}}$ to be the estimators at the last step, and $\mathcal{V}$ to be the matrix containing eigenvectors of $p^{-1}\bm{W}^{{\dagger}}(\bm{\Sigma}_{e}^{{\dagger}})^{-1}(\bm{W}^{{\dagger}})^{{\mathrm{\scriptscriptstyle T}}}$ corresponding to the descending eigenvalues. The maximum likelihood estimators are $\hat{\bm{W}}=\mathcal{V}^{{\mathrm{\scriptscriptstyle T}}}\bm{W}^{{\dagger}}$ and $\hat{\bm{\Sigma}}_{e}=\bm{\Sigma}_{e}^{{\dagger}}$. The estimator $\hat{\bm{S}}_{u}$ is obtained after estimating $\hat{\bm{W}}$ and $\hat{\bm{\Sigma}}_{e}$. As $\hat{\bm{S}}_{u}$ is not our focus in this method, we omit its derivation details in this paper.

With $\hat{\bm{W}}$ and $\hat{\bm{\Sigma}}_{e}$, we then estimate $\bm{U}_{i}$ using the generalized least squares estimator:

$$\hat{\bm{U}}_{i}=(\hat{\bm{W}}\hat{\bm{\Sigma}}_{e}^{-1}\hat{\bm{W}}^{\mathrm{\scriptscriptstyle T}})^{-1}\hat{\bm{W}}\hat{\bm{\Sigma}}_{e}^{-1}(\bm{X}_{i}-\bar{\bm{X}}).$$

(4)

Next, we treat the estimators $\hat{\bm{U}}_{1},\ldots,\hat{\bm{U}}_{n}$ as surrogate variables for the underlying unmeasured confounders to fit a high-dimensional generalized linear model. We then construct a debiased estimator for the parameter of interest by generalizing the decorrelated score method proposed by Ning and Liu (2017). Other debiasing approaches such as that of van de Geer et al. (2014) could also be similarly developed.

Initial $\ell_{1}$-Penalized Estimator: Recall that $\bm{\eta}=(\theta,\bm{v}^{{\mathrm{\scriptscriptstyle T}}},\bm{\beta}^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$. For vector $\bm{r}=(r_{1},\dots,r_{l})^{{\mathrm{\scriptscriptstyle T}}}$, define $\|\bm{r}\|_{1}=\sum_{j=1}^{l}|r_{j}|$. To fit a high-dimensional generalized linear model, we solve the $\ell_{1}$-penalized optimization problem

$$\hat{\bm{\eta}}=\underset{\bm{\eta}\in\mathbb{R}^{(p+K)}}{\operatorname{argmin}}\left[-\frac{1}{n}\sum_{i=1}^{n}\{y_{i}(\theta D_{i}+\bm{v}^{\mathrm{\scriptscriptstyle T}}\bm{Q}_{i}+\bm{\beta}^{\mathrm{\scriptscriptstyle T}}\hat{\bm{U}}_{i})-b(\theta D_{i}+\bm{v}^{\mathrm{\scriptscriptstyle T}}\bm{Q}_{i}+\bm{\beta}^{\mathrm{\scriptscriptstyle T}}\hat{\bm{U}}_{i})\}+\lambda(|\theta|+\|\bm{v}\|_{1})\right],$$

where $\hat{\bm{\eta}}=(\hat{\theta},\hat{\bm{v}}^{{\mathrm{\scriptscriptstyle T}}},\hat{\bm{\beta}}^{{\mathrm{\scriptscriptstyle T}}})^{\mathrm{\scriptscriptstyle T}}$ and $b(t)$ is a known function, given a specific generalized linear model. Throughout the manuscript, for notational convenience, let ${\bm{\zeta}}=(\bm{v}^{\mathrm{\scriptscriptstyle T}},\bm{\beta}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}$ be regression coefficients for the nuisance covariates and unmeasured confounders. Accordingly, we have $\bm{\eta}=(\theta,\bm{\zeta}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}$, with the goal of performing statistical inference on $\theta$. In the following, we construct a debiased estimator for $\theta$, generalizing the approach in Ning and Liu (2017) to situations involving unmeasured confounders.

A Debiased Estimator: Before we unfold the details of constructing a debiased estimator, we start with introducing some notations. Let $\bm{I}=E\{b^{\prime\prime}(\theta D_{i}+\bm{v}^{\mathrm{\scriptscriptstyle T}}\bm{Q}_{i}+\bm{\beta}^{\mathrm{\scriptscriptstyle T}}\bm{U}_{i})\bm{Z}_{i}\bm{Z}_{i}^{\mathrm{\scriptscriptstyle T}}\}$ be the Fisher information matrix, and let $\bm{w}^{\mathrm{\scriptscriptstyle T}}=\bm{I}_{\theta\bm{\zeta}}\bm{I}_{\bm{\zeta}\bm{\zeta}}^{-1}$, where $\bm{I}_{\theta\bm{\zeta}}$ and $\bm{I}_{\bm{\zeta}\bm{\zeta}}$ are corresponding block matrices of $\bm{I}$. In addition, let ${I}_{\theta\mid\bm{\zeta}}=E\{b^{\prime\prime}(\theta D_{i}+\bm{v}^{\mathrm{\scriptscriptstyle T}}\bm{Q}_{i}+\bm{\beta}^{\mathrm{\scriptscriptstyle T}}\bm{U}_{i})D_{i}(D_{i}-\bm{w}^{\mathrm{\scriptscriptstyle T}}\bm{M}_{i})\}$ be the partial Fisher information matrix, where $\bm{M}_{i}=(\bm{Q}_{i}^{\mathrm{\scriptscriptstyle T}},\bm{U}_{i}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}$ is a vector of nuisance covariates and unmeasured confounders. Finally, let $l(\theta,\bm{\zeta})=-n^{-1}\sum_{i=1}^{n}\{y_{i}(\theta D_{i}+\bm{\zeta}^{T}\dot{\bm{M}}_{i})-b(\theta D_{i}+\bm{\zeta}^{T}\dot{\bm{M}}_{i})\}$ with $\dot{\bm{M}}_{i}=(\bm{Q}_{i}^{\mathrm{\scriptscriptstyle T}},\hat{\bm{U}}_{i}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}$ be the loss function.

We define $S(\theta,\bm{\zeta})=\nabla_{\theta}l(\theta,\bm{\zeta})-\bm{w}^{\mathrm{\scriptscriptstyle T}}\nabla_{\bm{\zeta}}l(\theta,\bm{\zeta})$ as the generalized decorrelated score function, where $\nabla_{\theta}l(\theta,\bm{\zeta})$ and $\nabla_{\bm{\zeta}}l(\theta,\bm{\zeta})$ are the partial derivatives of the loss function with respect to $\theta$ and $\bm{\zeta}$, respectively. Different from the existing definition of the decorrelated score function in Ning and Liu (2017), the generalized decorrelated score function takes into account the effects induced by the unmeasured confounders. Specifically, in the presence of unmeasured confounders, the generalized decorrelated score function is uncorrelated with the score function corresponding to the nuisance covariates as well as the unmeasured confounders, i.e., $E\{S(\theta,\bm{\zeta})^{\mathrm{\scriptscriptstyle T}}\nabla_{\bm{\zeta}}l(\theta,\bm{\zeta})\}=\bm{I}_{\theta\bm{\zeta}}-\bm{I}_{\theta\bm{\zeta}}\bm{I}_{\bm{\zeta}\bm{\zeta}}^{-1}\bm{I}_{\bm{\zeta}\bm{\zeta}}=\bm{0}$. The debiased estimator of $\theta$ is constructed by solving for $\tilde{\theta}$ from the first-order approximation of the generalized decorrelated score function $\hat{S}(\hat{\theta},\hat{\bm{\zeta}})+\hat{I}_{\theta\mid\bm{\zeta}}(\tilde{\theta}-\hat{\theta})=0$. From the first-order approximation equation, we see that to establish the debiased estimator $\tilde{\theta}$, we need to construct two estimators $\hat{S}(\hat{\theta},\hat{\bm{\zeta}})$ and $\hat{I}_{\theta\mid\bm{\zeta}}$, and the key is to estimate $\bm{w}$.

We estimate $\bm{w}$ by solving the following convex optimization problem:

$$\hat{\bm{w}}=\underset{\bm{w}\in\mathbb{R}^{(p+K-1)}}{\operatorname{argmin}}~\frac{1}{2n}\sum_{i=1}^{n}\{\bm{w}^{\mathrm{\scriptscriptstyle T}}\nabla_{\bm{\zeta}\bm{\zeta}}l_{i}(\hat{\theta},\hat{\bm{\zeta}})\bm{w}-2\bm{w}^{{\mathrm{\scriptscriptstyle T}}}\nabla_{\bm{\zeta}\theta}l_{i}(\hat{\theta},\hat{\bm{\zeta}})\}+\lambda^{\prime}\|\bm{w}\|_{1},$$

where $l_{i}(\theta,\bm{\zeta})=-y_{i}(\theta D_{i}+\bm{\zeta}^{T}\dot{\bm{M}}_{i})+b(\theta D_{i}+\bm{\zeta}^{T}\dot{\bm{M}}_{i})$ is the $i$th component of the loss function and $\lambda^{\prime}>0$ is a sparsity tuning parameter for $\bm{w}$. Equivalently, the estimator $\hat{\bm{w}}$ is obtained by

$$\hat{\bm{w}}=\underset{\bm{w}\in\mathbb{R}^{(p+K-1)}}{\operatorname{argmin}}~\frac{1}{2n}\sum_{i=1}^{n}b^{\prime\prime}(\hat{\theta}D_{i}+\hat{\bm{\zeta}}^{{\mathrm{\scriptscriptstyle T}}}\dot{\bm{M}}_{i})(D_{i}-\bm{w}^{{\mathrm{\scriptscriptstyle T}}}\dot{\bm{M}}_{i})^{2}+\lambda^{\prime}\|\bm{w}\|_{1}.$$

(5)

The estimator $\hat{\bm{w}}$ is constructed with the intuition of finding a sparse vector such that the generalized decorrelated score function is approximately zero. This coincides with the intuition to solve for $\tilde{\theta}$ from the first-order approximation of the generalized decorrelated score function. Under the null hypothesis $H_{0}:\theta^{*}=\theta^{0}$, we estimate the generalized decorrelated score function and the partial Fisher information matrix by

	$\displaystyle\hat{S}(\theta^{0},\hat{\bm{\zeta}})$	$\displaystyle=$	$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}\{y_{i}-b^{\prime}(\theta^{0}D_{i}+\hat{\bm{v}}^{\mathrm{\scriptscriptstyle T}}\bm{Q}_{i}+\hat{\bm{\beta}}^{\mathrm{\scriptscriptstyle T}}\hat{\bm{U}}_{i})\}(D_{i}-\hat{\bm{w}}^{\mathrm{\scriptscriptstyle T}}\dot{\bm{M}}_{i});$
	$\displaystyle\hat{{I}}_{\theta\mid\bm{\zeta}}$	$\displaystyle=$	$\displaystyle\frac{1}{n}\sum_{i=1}^{n}b^{\prime\prime}(\hat{\theta}D_{i}+\hat{\bm{v}}^{\mathrm{\scriptscriptstyle T}}\bm{Q}_{i}+\hat{\bm{\beta}}^{\mathrm{\scriptscriptstyle T}}\hat{\bm{U}}_{i})D_{i}(D_{i}-\hat{\bm{w}}^{\mathrm{\scriptscriptstyle T}}\dot{\bm{M}_{i}}).$

The debiased estimator can then be constructed as $\tilde{\theta}=\hat{\theta}-(\hat{I}_{\theta\mid\bm{\zeta}})^{-1}\hat{S}(\theta^{0},\hat{\bm{\zeta}})$. We will show in Section 4 that the debiased estimator $\tilde{\theta}$ is asymptotically normal. Subsequently, the $(1-\alpha)\times 100\%$ confidence interval for $\theta$ can be constructed as

$$\{\tilde{\theta}-(n\hat{I}_{\theta\mid\bm{\zeta}})^{-1/2}{\Phi}^{-1}(1-{\alpha}/2),\ \tilde{\theta}+(n\hat{I}_{\theta\mid\bm{\zeta}})^{-1/2}{\Phi}^{-1}(1-{\alpha}/2)\},$$

(6)

where $\Phi(t)$ is the cumulative distribution function for the standard normal random variable.

Recall that our proposed method yields estimators $\hat{\bm{\eta}}$, $\hat{\bm{w}}$, and the debiased estimator $\tilde{\theta}$. In this section, we first establish upper bounds for the estimation errors of $\hat{\bm{\eta}}$ and $\hat{\bm{w}}$ under the $\ell_{1}$ norm. Subsequently, we show that the debiased estimator $\tilde{\theta}$ is asymptotically normal.

For a vector $\bm{r}=(r_{1},\dots,r_{l})^{{\mathrm{\scriptscriptstyle T}}}$, let $\|\bm{r}\|_{q}=(\sum_{j=1}^{l}|r_{j}|^{q})^{1/q}$ for $q\geq 1$ and let $\|\bm{r}\|_{\infty}=\max_{j=1,\ldots,l}|r_{j}|$. For any matrix $\bm{A}$, let $\lambda_{\max}(\bm{A})$ and $\lambda_{\min}(\bm{A})$ represent the largest and smallest eigenvalues of $\bm{A}$, respectively. Moreover, for sequences $\{a_{n}\}$ and $\{b_{n}\}$, we write $a_{n}\lesssim b_{n}$ if there exists a constant $C>0$ such that $a_{n}\leq Cb_{n}$ for all $n$, and $a_{n}\asymp b_{n}$ if $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$. For a sub-exponential random variable $Y_{1}$, we write $\|Y_{1}\|_{\varphi_{1}}=\inf[s>0:E\{\exp(Y_{1}/s)\}\leq 2]$ as the sub-exponential norm. For a sub-Gaussian random variable $Y_{2}$, we write $\|Y_{2}\|_{\varphi_{2}}=\inf[s>0:E\{\exp{(Y_{2}^{2}/s^{2})}\}\leq 2]$ as the sub-Gaussian norm. Throughout the manuscript, we will use an asterisk on the upper subscript to indicate the population parameters. In addition, we define $s_{\eta}=\text{card}\{(j:\eta^{*}_{j}\neq 0)\}$ and $s_{w}=\text{card}\{(j:\bm{w}^{*}_{j}\neq 0)\}$ as the cardinalities of $\bm{\eta}^{*}$ and $\bm{w}^{*}$, respectively. All of our theoretical analysis are performed under the regime in which $n$, $p$, $s_{\eta}$, and $s_{w}$ are allowed to increase, and the number of unmeasured confounders $K$ is fixed.

We start with some conditions on the factor model in (2). Similar conditions were also considered in Bai and Li (2016) in the context of high-dimensional approximate factor model.

Assumption 1 is more general than assumptions in classical factor analysis (Anderson and Amemiya, 1988; Bai, 2003; Fan et al., 2013; Bai and Li, 2016). Instead of constraining all the $\bm{E}_{i}$ to have a diagonal covariance matrix, we now only require the higher-order moment of $\bm{E}_{i}$ to be bounded, the diagonal entries $\sigma_{j}^{*}$’s to be bounded, as well as the magnitudes of the correlations among entries of $\bm{E}_{i}$ to be controlled in Assumption 1. The conditions are mild as we only control the magnitude of correlations rather than assuming zero correlations as in classical factor analysis. Assumption 1(e) is a regularity condition for restricting the parameters. Overall, Assumption 1 follows standard conditions for the approximate factor model in Bai and Li (2016) and is required for the estimation consistency of the unmeasured confounders.

Comparing Assumption 1 to the dense confounding assumption (A2) imposed on the linear framework in Guo et al. (2022), our assumption is mild as it holds for a broad regime of $n$ and $p$. Specifically, the dense confounding assumption (A2) in Guo et al. (2022) is related to our assumption 1(e) that $\lim_{p\rightarrow\infty}p^{-1}\bm{W}^{*}(\bm{\Sigma}_{e}^{*})^{-1}(\bm{W}^{*})^{\mathrm{\scriptscriptstyle T}}=\bm{\Gamma}^{*}$ where $\bm{\Gamma}^{*}$ is positive definite. Our assumption follows from classical factor analysis literature and as a result, we have $\lambda_{q}(\bm{W}^{*})\asymp\sqrt{p}$, where $\lambda_{q}(\bm{W}^{*})$ is the $q$-th singular value of the factor loading matrix $\bm{W}^{*}$. With $\bm{\gamma}^{{\mathrm{\scriptscriptstyle T}}}=(\theta,\bm{v}^{{\mathrm{\scriptscriptstyle T}}})$ being the coefficient for all covariates and $\gamma_{j}$ being the individual coefficient of interest, the dense confounding assumption mainly requires that

$$\lambda_{q}(\bm{W}_{-j}^{*})\gg\max\left\{M\sqrt{Kpn^{-1}}(\log p)^{3/4},\sqrt{MK}p^{1/4}(\log p)^{3/8},\sqrt{Kn\log p}\right\},$$

where $\bm{W}_{-j}^{*}$ denotes the factor loading matrix $\bm{W}$ with $j$-th column removed. Guo et al. (2022) focuses on the setting $p\gg n$ and they point out that in the high-dimensional regime and under certain settings, the dense confounding assumption holds with high probability. Specifically, when the entries of $\bm{W}$ are i.i.d. Sub-Gaussian with zero mean and variance $\sigma_{w}^{2}$, it holds that $\lambda_{q}(\bm{W}_{-j})\asymp\sqrt{p}\sigma_{w}$ and the dense confounding assumption requires $p\gg Kn\log p$ and $\min\{n,p\}\gg K^{3}(\log p)^{3/2}M^{2}$ to make $\sigma_{w}$ diminish to zero. However, this condition is restricted to the high-dimensional regime and may not hold when $p$ is of relatively lower order.

Moreover, without loss of generality, we assume a working identifiability condition: $\bm{S}_{u}=\bm{I}_{K}$ and $p^{-1}\bm{W}^{*}(\bm{\Sigma}_{e}^{*})^{-1}(\bm{W}^{*})^{\mathrm{\scriptscriptstyle T}}$ is a diagonal matrix with distinct entries. Note that the aforementioned working identifiability condition is for presentation purpose only and not an assumption on the model structure. As presented in Appendix B, when such a working identifiability condition is not satisfied, our theoretical results in Theorems 2 and 3 are still valid. We also illustrate this via simulation in Section 5.1. Next, we impose some assumptions on the generalized linear model with unmeasured confounders in (1).

In the absence of unmeasured confounders, similar conditions in Assumption 2 can be implied from the conditions in Theorem 3.3 in van de Geer et al. (2014) and Assumption E.1 in Ning and Liu (2017). When $\bm{U}_{i}$ and $\bm{X}_{i}$ are binary or categorical, Assumption 2(b) holds with a constant $M>0$. When $\bm{U}_{i}$ and $\bm{X}_{i}$ are sub-exponential random vectors, Assumption 2(b) holds with $M=cn^{-1/2}(\log p)^{1/2}$ for some constant $c>0$, with probability at least $1-p^{-1}$. Assumption 2(d) imposes mild regularity conditions on the function $b(t)$, and is commonly used in analyzing high-dimensional generalized linear models without unmeasured confounders. We require the function $b(t)$ to be at least twice differentiable and $b^{\prime}(t)$ and $b^{\prime\prime}(t)$ to be bounded. Specifically, $|b^{\prime\prime}(t_{1})-b^{\prime\prime}(t)|\leq B|t_{1}-t|b^{\prime\prime}(t)$ can be implied by $|b^{\prime\prime\prime}(t)|\leq Bb^{\prime\prime}(t)$ when $b^{\prime\prime\prime}(t)$ exists, which is a weaker condition than the self-concordance property (Bach, 2010). This boundary assumption is important for the concentration of the Hessian matrix of the loss function. Assumption 2(d) holds for commonly used generalized linear models. For example, for logistic regression where $b(t)=\log\{1+\exp(t)\}$, Assumption 2(d) holds with $B=1$, and $a_{1},a_{2}$ extended to infinity and it can be similarly verified to hold at $B=\max\{1,\exp(a_{2}+\epsilon)\}$ for Poisson regression and $B=\max\{1/(a_{2}+\epsilon)^{2},2/|a_{2}+\epsilon|\}$ for exponential regression with $a_{2}<0$. For linear model, Assumption 2 can be relaxed as stated in Remark 4.

The theoretical analysis on the unmeasured confounders estimator is important in establishing the theoretical guarantee for our debiased method. As the decomposition techniques commonly used in linear models may not be applicable in generalized linear model settings, it is necessary to establish more general and refined intermediate results as the foundation of our theoretical analysis (see Remark 8 in Appendix D.1 for more details). We first present a uniform convergence result for the estimators of the unmeasured confounders.

From Proposition 1, the estimator of unmeasured confounders $\hat{\bm{U}}_{i}$ uniformly converges to $\bm{U}_{i}$ at $p\gg\log n$, which holds naturally under the high-dimensional regime $p\gg n$. Moreover, the convergence rate $O_{p}(n^{-1/2}+p^{-1/2}$ $(\log n)^{1/2})$ is of a similar order to the convergence rates of principal component estimators (Fan et al., 2013; Wang and Fan, 2017). The estimation results of $\hat{\bm{\eta}}$ and $\hat{\bm{w}}$ are dependent on the accuracy of the estimator of unmeasured confounders and we next establish upper bounds on the estimation errors for $\hat{\bm{\eta}}$ and $\hat{\bm{w}}$.

The effect of unmeasured confounders enters the rate of the estimation error in our results through the term $p^{-1/2}(\log n)^{1/2}$. Under the high-dimensional setting with $p\gg n$, the unmeasured confounders effect is dominated by $n^{-1/2}(\log p)^{1/2}$ and as a result, the convergence rates in Theorem 2 are of the same order as in the oracle case when the unmeasured confounders are assumed to be known (van de Geer et al., 2014; Ning and Liu, 2017). With the estimation consistency established, we show that the debiased estimator from the proposed estimation method is asymptotically normal, and thus valid confidence intervals can be constructed.

The result in Theorem 3 is similar to existing inference results for high-dimensional generalized linear models without unmeasured confounders (van de Geer et al., 2014; Ning and Liu, 2017; Ma et al., 2021; Shi et al., 2021; Cai et al., 2023). Theorem 3 is established under the condition that the estimation errors of $\hat{\bm{\eta}}$ and $\hat{\bm{w}}$ in Theorem 2 are $o_{p}(1)$. As a consequence of the asymptotic normality result in Theorem 3 and that $\hat{{I}}_{\theta\mid\bm{\zeta}}$ is consistent estimator for ${I}_{\theta\mid\bm{\zeta}}^{*}$, the construction of the confidence interval in (6) is valid.

To illustrate the performance of our proposed method, we conduct numerical experiments including simulation studies and an application of our method to a genetic data set.