G-HIVE: Parameter Estimation and Approximate Inference for Multivariate Response Generalized Linear Models with Hidden Variables

TL;DR

This paper introduces G-HIVE, a unified framework for parameter estimation and approximate inference in multivariate response GLMs with hidden variables, addressing challenges in non-linear models.

Contribution

It proposes a novel quasi-likelihood method with bias correction and an approximate inference framework for models with hidden variables, extending existing linear regression approaches.

Findings

G-HIVE effectively estimates parameters in complex models.

The method reduces bias and improves inference accuracy.

Simulations and real data show G-HIVE outperforms baseline methods.

Abstract

In practice, there often exist unobserved variables, also termed hidden variables, associated with both the response and covariates. Existing works in the literature mostly focus on linear regression with hidden variables. However, when the regression model is non-linear, the presence of hidden variables leads to new challenges in parameter identification, estimation, and statistical inference. This paper studies multivariate response generalized linear models (GLMs) with hidden variables. We propose a unified framework for parameter estimation and statistical inference called G-HIVE, short for 'G'eneralized - 'HI'dden 'V'ariable adjusted 'E'stimation. Specifically, based on factor model assumptions, we propose a modified quasi-likelihood approach to estimate an intermediate parameter, defined through a set of reweighted estimating equations. The key of our approach is to construct the proper weight, so that the first-order asymptotic bias of the estimator can be removed by orthogonal projection. Moreover, we propose an approximate inference framework for uncertainty quantification. Theoretically, we establish the first-order and second-order asymptotic bias and the convergence rate of our estimator. In addition, we characterize the accuracy of the Gaussian approximation of our estimator via the Berry-Esseen bound, which justifies the validity of the proposed approximate inference approach. Extensive simulations and real data analysis results show that G-HIVE is feasibly implementable and can outperform the baseline method that ignores hidden variables.

Figures

Tables

Table 1. Table 1 : Inference results showing the length of the confidence interval, standard error of the estimates, and the coverage probabilities for u ⊤ ​ P B ⟂ ​ F ∗ ​ v u^{\top}P^{\perp}_{B}F^{*}v and u ⊤ ​ Θ ​ v u^{\top}\Theta v at level α = 0.05 \alpha=0.05 for u ⊤ = v ⊤ = [ 1 , 0 , 0 , 0 ] u^{\top}=v^{\top}=[1,0,0,0] , averaged over r = 100 r=100 repetitions.

	$n$	CI Length	Std. Err.	Coverage for ${\left(P_B^{⟂}{F}^{\ast }\right)}_{11}$	Coverage for ${\Theta }_{11}$
data driven g-hive	40	7.77	1.98	0.98	0.98
	70	6.40	1.63	0.99	0.99
naive mle	40	1.75	0.45	0.77	0.75
	70	1.11	0.28	0.62	0.61

Table 2. Table 2 : The coefficient values relating the explanatory variables to the main response variable (diagnosis of mesothelioma) obtained with naive mle and G-hive in the lung cancer dataset.

Method	Lung Side	WCC	Platelets	Sedim.	Albumin	Glucose	PLD	Pleural Prot.	Pleural Thick.
naive mle	0.2057	-0.1091	-0.2896	0.0895	0.1009	0.0613	-0.0496	-0.1503	0.0861
G-hive	0.1941	-0.1116	-0.4111	0.1759	0.1148	0.0706	0.0889	-0.1395	0.1496

Table 3. Table 3 : G-hive results for the coefficient values relating the covariates to the response variables in the NHANES dataset with n = 320 n=320 observations, M = 5 M=5 binary response variables and p = 5 p=5 observed covariates.

	Age	Income	SysBP	BMI	Glucose
Depression	0.0240	-0.4243	-0.2453	0.5891	-0.7752
Diabetes	0.3589	0.1676	0.1695	0.4284	0.9332
Osteoporosis	-0.0595	0.0475	0.0286	-0.1754	0.0643
Hypertension	0.3733	0.1202	0.4070	0.5687	0.6886
Arthritis	0.3566	-0.0667	-0.2127	0.7989	-0.3615

Table 4. Table 4 : Results of naive mle and G-hive on the reduced model. The estimates in the second column represent the effect of “Income” on the response variables without taking into account the effect of “Age”.

Outcome	Intercept	Income
Smoke	$-1.45$	$-0.63$
Diabetes	$-1.60$	$-0.17$
Hypertension	$-0.65$	$0.16$
Arthritis	$-1.14$	$0.19$
Depression	$-1.31$	$-0.33$

Table 5. (a) Results of naive mle on the reduced model.

Outcome	Intercept	Income
Smoke	$-1.45$	$-0.63$
Diabetes	$-1.60$	$-0.17$
Hypertension	$-0.65$	$0.16$
Arthritis	$-1.14$	$0.19$
Depression	$-1.31$	$-0.33$

Table 6. (b) Results of G-hive on the reduced model.

Outcome	Intercept	Income
Smoke	$-1.90$	$-0.47$
Diabetes	$-1.02$	$-0.38$
Hypertension	$-0.45$	$0.09$
Arthritis	$-1.05$	$0.14$
Depression	$-1.23$	$-0.35$

Table 7. Table 5 : Results of naive mle on the full model. The estimates in the third column represent the effect of “Income” on the response variables while taking into account the effect of “Age”.

Outcome	Intercept	Age	Income
Smoke	$-1.52$	$-0.51$	$-0.54$
Diabetes	$-1.91$	$1.03$	$-0.36$
Hypertension	$-0.73$	$0.73$	$0.04$
Arthritis	$-1.32$	$0.87$	$0.07$
Depression	$-1.32$	$-0.23$	$-0.29$

Equations

$$\displaystyle f(Y_{m}|X,Z)=\exp\big{[}\{Y_{m}\cdot(\Theta_{m}X+B_{m}Z)-b(\Theta_{m}X+B_{m}Z)\}/\phi\leavevmode\nobreak\ +\leavevmode\nobreak\ c(Y_{m},\phi)\big{]}$$

$$X=AZ+W,$$

$$\displaystyle F_{m}^{\textsc{MLE}}\leavevmode\nobreak\ =\leavevmode\nobreak\ \underset{F_{m}\in\mathbb{R}^{p}}{\text{arg max}}\leavevmode\nobreak\ \mathbb{E}\Big{\{}Y_{m}\cdot(F_{m}X)-b(F_{m}X)\Big{\}}.$$

$$\displaystyle F_{m}^{\textsc{MLE}}-\Theta_{m}=\mathbb{E}\Big{\{}b^{\prime\prime}(\Theta_{m}X+B_{m}Z)B_{m}ZZ^{T}\Big{\}}A^{T}\Big{\{}\mathbb{E}(b^{\prime\prime}(\Theta_{m}X+B_{m}Z)XX^{T})\Big{\}}^{-1}+Rem_{m},$$

$$P_B^{\perp }{F}^{\text{textsc}MLE}\approx P_B^{\perp }(\Theta +\mathbb{E}(DBZ{Z}^T)A^T\mathbb{E}(DX{X}^T))\ne \Theta ,$$

$$\displaystyle\mathbb{E}\bigg{[}\bigg{\{}\frac{Y_{m}-b^{\prime}(F_{m}^{*}X)}{b^{\prime\prime}(F_{m}^{*}X)}\bigg{\}}X^{T}\bigg{]}=0.$$

$$\displaystyle F_{m}^{*}-\Theta_{m}=B_{m}\mathbb{E}(ZZ^{T})A^{T}\Big{\{}\mathbb{E}(XX^{T})\Big{\}}^{-1}+Rem_{m}^{\prime},$$

$$P_B^{\perp }{F}^{\ast }\approx P_B^{\perp }(\Theta +B_m\mathbb{E}(Z{Z}^T)A^T\{\mathbb{E}(X{X}^T){\}}^{-1})=\Theta .$$

$${\overline{ϵ}}_m:$$

$$\displaystyle\mathbb{E}(\bar{\epsilon}\bar{\epsilon}^{T})=\mathbb{E}({\epsilon}{\epsilon}^{T})+B\mathbb{E}\Big{[}(Z-\Gamma X)(Z-\Gamma X)^{T}\Big{]}B^{T}+Rem^{\prime\prime},$$

$$\displaystyle\epsilon_{m}:=\frac{Y_{m}-b^{\prime}(\Theta_{m}X+B_{m}Z)}{b^{\prime\prime}(\Theta_{m}X+B_{m}Z)},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \epsilon:=\big{[}{\epsilon}_{1},...,{\epsilon}_{M}\big{]}^{T}$$

$${\hat{F}}_m^{(D_1)}=\arg \max Q_m^{(D_1)}(F_m),\text{leavevmode}\text{leavevmode}\text{where}\text{leavevmode}\text{leavevmode}{Q}_m^{(D_1)}(F_m)=\frac{1}{|D_1|}\sum _{i\in D_1}{\int }_0^{F_m{X}^{(i)}}\frac{Y_m^{(i)}-b^{\prime }(\eta )}{b^{\prime \prime }(\eta )}d\eta .$$

$${\widehat{ϵ}}_m^{(i)}=\frac{Y_m^{(i)}-b^{\prime }({\hat{F}}_m^{(D_1)}{X}^{(i)})}{b^{\prime \prime }({\hat{F}}_m^{(D_1)}{X}^{(i)})}.$$

$${\hat{\Sigma }}^{(D_2)}=\frac{1}{|D_2|}\sum _{i\in D_2}{\widehat{ϵ}}^{(i)}({\widehat{ϵ}}^{(i)}{)}^T,$$

$${\hat{F}}_m=({\hat{F}}_m^{(D_1)}+{\hat{F}}_m^{(D_2)})/2,\text{leavevmode}\text{leavevmode}\text{and}\text{leavevmode}\text{leavevmode}\hat{\Sigma }=({\hat{\Sigma }}^{(D_1)}+{\hat{\Sigma }}^{(D_2)})/2.$$

$$\hat{K}\text{leavevmode}=\text{leavevmode}\underset{j\in \{1,2,...,\bar{K}\}}{\mathrm{argmax}}\text{leavevmode}\frac{{\hat{\lambda }}_j}{{\hat{\lambda }}_{j+1}}$$

$$\big{|}\big{|}F^{*}_{m}-\Theta_{m}\big{|}\big{|}_{2}=O\bigg{(}\frac{1}{\sqrt{p}}\bigg{)}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \big{|}\big{|}(P_{B}^{\perp}F^{*})_{m}-\Theta_{m}\big{|}\big{|}_{2}=O\bigg{(}\frac{1}{p}\bigg{)},$$

$$\sqrt{n}({\hat{F}}_{mj}-F_{mj}^{\ast })=\frac{1}{\sqrt{n}}\sum _{i=1}^n\frac{Y_m^{(i)}-b^{\prime }(F_m^{\ast }{X}^{(i)})}{b^{\prime \prime }(F_m^{\ast }{X}^{(i)})}{e}_j^T{G}_m^{-1}{X}^{(i)}+o_p(1),$$

$$\sqrt{n}(\hat{F}_{mj}-\Theta_{mj})=\sqrt{n}(\hat{F}_{mj}-F^{*}_{mj})+\sqrt{n}(F^{*}_{mj}-\Theta_{mj})\rightarrow_{d}N(0,\sigma^{2}_{mj})+O_{p}\Big{(}\sqrt{\frac{n}{p}}\Big{)},$$

$$G_{m}=\mathbb{E}\Big{(}1+\zeta^{(i)}_{m}(F^{*}_{m})\Big{)}X^{(i)}X^{(i)T},$$

$${\zeta }_m^{(i)}(F_m)=\frac{(Y_m^{(i)}-b^{\prime }(F_m{X}^{(i)}))b^{\prime \prime \prime }(F_m{X}^{(i)})}{\{b^{\prime \prime }(F_m{X}^{(i)}){\}}^2}.$$

$$\sqrt{n}{u}^T(\hat{\Theta }-P_B^{\perp }{F}^{\ast })v/(s_n/\sqrt{n}){\to }_{d}N(0,1),$$

$$\hat{G}_{m}=\frac{1}{n}\sum_{i=1}^{n}\Big{(}1+\zeta^{(i)}_{m}(\hat{F}_{m})\Big{)}X^{(i)}X^{(i)T}$$

$$(u^T\hat{\Theta }v-q_{1-\alpha /2}{\hat{s}}_n/\sqrt{n},u^T\hat{\Theta }v+q_{1-\alpha /2}{\hat{s}}_n/\sqrt{n}).$$

$$F_m^{\ast }-{\Theta }_m=B_m{\Sigma }_Z{A}^T{\Sigma }_X^{-1}+Re{m}_m^{\prime },$$

$$\displaystyle\max_{1\leq m\leq M}\mathbb{E}\Big{[}\big{(}\Theta_{m}X+B_{m}Z-F_{m}^{*}X\big{)}^{4}\Big{]}=O\bigg{(}\frac{1}{p^{2}}\bigg{)}$$

$$\displaystyle\max_{1\leq m\leq M}||F_{m}^{*}-\Theta_{m}||_{2}=O\bigg{(}\frac{1}{\sqrt{p}}\bigg{)}.$$

$$\displaystyle\frac{1}{\sqrt{M}}\big{|}\big{|}F^{*}-\Theta\big{|}\big{|}_{F}=O\bigg{(}\frac{1}{\sqrt{p}}\bigg{)}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \text{{and}}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \frac{1}{\sqrt{M}}\big{|}\big{|}P_{B}^{\perp}F^{*}-\Theta\big{|}\big{|}_{F}=O\bigg{(}\frac{1}{p}\bigg{)}.$$

$$\displaystyle\max_{1\leq m\leq M}\big{|}\big{|}\hat{F}_{m}-F_{m}^{*}\big{|}\big{|}_{2}=O_{p}\Bigg{(}\sqrt{\frac{p\log(M\vee p)}{n}}\Bigg{)}.$$

$$({\hat{F}}_m-F_m^{\ast })v=\frac{1}{n}\sum _{i=1}^n\frac{Y_m^{(i)}-b^{\prime }(F_m^{\ast }{X}^{(i)})}{b^{\prime \prime }(F_m^{\ast }{X}^{(i)})}{v}^T{G}_m^{-1}{X}^{(i)}+Re{m}_m^{\prime \prime },$$

Full text

Keywords: Generalized linear models, multivariate response data, non-linear regression, hidden variables, unmeasured confounders, parameter estimation, approximate inference.

1 Introduction

In many regression problems, due to measurement limitations or ethical considerations, there often exist unobserved variables, also referred to as hidden variables. For example, in the analysis of high-throughput genomic data, researchers have long been aware that the measurements can be affected by many unobserved factors such as laboratory conditions, preparation time, and reagent lots (Irizarry et al., 2005; Luo and Wei, 2019). These factors are called batch effects, which can be modeled as hidden variables (Leek and Storey, 2007). Similarly, in biomedical studies, the onset of a disease is likely associated with several unmeasured variables, such as environmental factors or habitual patterns (Katsaouni et al., 2021). Ignoring hidden variables in the statistical analysis may introduce estimation bias and potentially lead to misleading scientific conclusions. Therefore, there is a pressing need to develop statistical methods that deal with hidden variables in a general regression framework, and that in particular, are applicable to binary or categorical data.

This paper studies the multivariate response generalized linear model with hidden variables. Specifically, we assume that the $M$-dimensional response variable $Y=(Y_{1},...,Y_{M})^{T}$ given the observed covariates $X\in\mathbb{R}^{p}$ and hidden variables $Z\in\mathbb{R}^{K}$ follows the generalized linear model (GLM) with the canonical link

[TABLE]

where $b(\cdot)$ and $c(\cdot)$ are known functions and $\phi$ is the dispersion parameter. The parameters $\Theta_{m}$ and $B_{m}$ are the $m$-th row of coefficient matrices $\Theta\in\mathbb{R}^{M\times p}$ and $B\in\mathbb{R}^{M\times K}$, respectively. Given $n$ i.i.d copies of $(Y,X)$, we are interested in the estimation and inference of the coefficient matrix $\Theta$, the association between $X$ and $Y$ in the presence of hidden variables $Z$. In this work, we consider the regime where $p,M,K$ are all allowed to grow with the sample size $n$, where $p\leq n$ and $K\leq M$ hold.

In this work, we propose a unified framework for parameter estimation and statistical inference called G-hive, short for Generalized - HIdden Variable adjusted Estimation. Since the hidden variable $Z$ is random and unobserved, the coefficient matrix $\Theta$ is generally not identifiable. To make $\Theta$ (asymptotically) identifiable and estimable, we impose a factor model in (2) that relates $X$ and $Z$ (Bai, 2003; Fan et al., 2013, 2008). However, under these model assumptions, the distribution of $Y$ given $X$ does not follow a GLM and is indeed intractable, since we do not impose any parametric assumption on the distribution of $Z$. To overcome this challenge, we carefully construct a modified quasi-likelihood for a new estimand $F^{*}$, which is defined through a set of reweighted estimating equations. The rationale behind the reweighted estimating equations is that the resulting first-order approximation of the bias $F^{*}-\Theta$ is shown to belong to the column space of $B$, which can be removed by estimating the projection matrix $P_{B}=B(B^{T}B)^{-1}B^{T}$. The intuition of our approach is explained in Section 2.2. Under certain identifiability conditions, using $F^{*}$ as a bridge, we introduce the estimator $\hat{\Theta}=\hat{P}_{B}^{\perp}\hat{F}$, where $\hat{P}_{B}$ is obtained by applying PCA to a carefully constructed weighted covariance matrix and $\hat{F}$ is the maximum modified quasi-likelihood estimator. Since our approach yields a tight pipeline, it is straightforward to implement in practice.

In addition, we propose an approximate inference framework for uncertainty quantification. Unlike classical inference results for models without hidden variables, a new, unpleasant phenomenon of the estimator $\hat{\Theta}$ is that the asymptotic bias may dominate the stochastic error. Consequently, the limiting distribution of the estimator $\hat{\Theta}$ is no longer centered at $\Theta$ with the $\sqrt{n}-$rate. To address this issue, we shift the target parameter from $\Theta$ to $P_{B}^{\perp}F^{*}$ (or $F^{*}$), which corresponds to the second-order (or first-order) approximation of $\Theta$. Intuitively, $P_{B}^{\perp}F^{*}$ can be viewed as the correct limiting value of $\hat{\Theta}$, and therefore a confidence interval based on the limiting distribution of $\hat{\Theta}$ yields the desired coverage probability for $P_{B}^{\perp}F^{*}$. For this reason, we refer to this approach as second-order approximate inference.

Theoretically, our first key result shows that the approximation bias satisfies $\|F^{*}-\Theta\|_{F}/\sqrt{M}=O(1/\sqrt{p})$ and $\|P_{B}^{\perp}F^{*}-\Theta\|_{F}/\sqrt{M}=O(1/p)$. An interesting implication of this result is that collecting more observed covariates can mitigate the approximation bias. Moreover, it also explains why $P_{B}^{\perp}F^{*}$ (or $F^{*}$) is called the second-order (or first-order) approximation of $\Theta$ in our inference framework. Next, we establish the convergence rate of the estimator $\hat{\Theta}$. In particular, we show that, under mild conditions (e.g., $M$ is large enough), the convergence rate of $\hat{\Theta}$ is faster than that of $\hat{F}$, which corresponds to the baseline naive estimator that ignores the hidden variables. Finally, we characterize the accuracy of the Gaussian approximation of $\hat{\Theta}$ via the Berry–Esseen bound, which justifies the validity of the proposed approximate inference approach.

1.1 Related Literature

This work is most related to surrogate variable analysis (SVA) proposed by Leek and Storey (2007), and more recently developed by Lee et al. (2017); Wang et al. (2017); McKennan and Nicolae (2019); Bing et al. (2022, 2023), the last of which proposed a novel factor model based bias correction approach for multivariate response linear regression with hidden variables, which is a special case of GLMs. However, their approach is only applicable to linear regression. The challenge of extending their approach to GLMs is detailed in Section 2.1. Compared to these works, our main methodological novelty is that we propose to calibrate the residual by an approximate inverse variance weighting scheme. Such a calibration step is essential under the GLM for parameter identification, estimation consistency, and asymptotic normality. Theoretically, we discover a unique result in that our estimator under the GLM inherently has an asymptotic bias which decreases with $p$ but may still dominate the stochastic error. As a result, there is an interesting and much more delicate interplay between $p$ and $M$ in both the estimation error and the Berry–Esseen bound for Gaussian approximation.

Along this line, a recent work by Du et al. (2025) studied simultaneous inference with unmeasured confounders when $p\gg n$. While they focused on the same GLM as in our (1), their imposed model for the unmeasured confounders is different from our factor model in (2), and consequently the corresponding assumptions on their model are different from our Assumption 3. In particular, our theory is established under a more challenging setting, where the covariance matrix of $X$ has spiked eigenvalues. Their proposed method is a joint maximum likelihood approach, which requires the estimation of all coefficient matrices as well as the latent factors $Z$ for each sample. In contrast, our method is computationally more convenient and avoids estimating the unknown factors.

Another recent direction of relevance, grouped together under the term spectral deconfounding, includes work by Ćevid et al. (2020); Guo et al. (2022); Fan et al. (2024); Wang and Shah (2025); Sun et al. (2024) among others, and considers estimation and inference in high-dimensional regression models with unmeasured confounders. For example, Ouyang et al. (2023) focused on inference in high-dimensional GLMs with unmeasured confounders by generalizing the decorrelated score approach from Ning and Liu (2017) to account for the effects induced by the unmeasured confounders. This work is similar to ours, but fundamentally different in that their response is assumed to be univariate. In contrast, we show that with multiple response variables, we can estimate the parameters in a collaborative way, improving the convergence rate compared to the univariate case.

Alternatively, one may view hidden variables as random effects or latent factors. The usage of random effects or latent variables in GLMs have many different forms in the literature (Bartholomew et al., 2011; McCulloch, 2001). For example, Huber et al. (2004) introduced generalized linear latent variable models without any observed covariates, and proposed a Laplacian approximation to estimate the coefficient of the latent variable. A similar approach was also considered for generalized linear mixed effect models (Breslow and Clayton, 1993). All these works differ substantially from our approach.

1.2 Notation

For any vector $v\in\mathbb{R}^{d}$ and some real number $q\geq 0$, we define its $L_{q}$ norm as $||v||_{q}=(\sum_{j=1}^{d}|v_{j}|^{q})^{1/q}$. For any matrix $H\in\mathbb{R}^{d_{1}\times d_{2}}$, we denote by $||H||_{\text{op}}$ and $||H||_{F}$ the operator norm and the Frobenius norm, respectively. $||H||_{\infty}=\max_{i}\sum_{j}|h_{ij}|$ denotes the maximum absolute row sum. Following the notation in Vershynin (2018), for any sub-Gaussian random variable (or vector) $h_{2}$, let $||h_{2}||_{\psi_{2}}$ denote its sub-Gaussian norm, and for any sub-exponential random variable (or vector) $h_{1}$, let $||h_{1}||_{\psi_{1}}$ denote its sub-exponential norm. For any symmetric matrix $H$, we write $\lambda_{k}(H)$ to denote its $k$-th largest eigenvalue, and $\lambda_{\min}(H)$ and $\lambda_{\max}(H)$ for the smallest and largest eigenvalues, respectively. For any two sequences $a_{n}$ and $b_{n}$, we write $a_{n}\lesssim b_{n}$ if there exists some fixed positive constant $C$ such that $a_{n}\leq Cb_{n}$. We also use the following notation to refer to the maximum and minimum: $a\vee b=\max(a,b)$, $a\wedge b=\min(a,b)$.

2 Informal Analysis of Parameter Identifiability

In this section, we first introduce our model setup and highlight the challenges of model identifiability under the GLM compared to the linear case, and afterwards we offer the intuition of our proposed approach.

2.1 Model Setup and Background

Recall that given the observed covariates $X\in\mathbb{R}^{p}$ and hidden variables $Z\in\mathbb{R}^{K}$, the response variable $Y=(Y_{1},...,Y_{M})^{T}$ follows the GLM (1). For simplicity, we assume that $\phi$ is known (set $\phi=1$) and that $X$ and $Z$ have zero mean. The parameters $\Theta_{m}$ and $B_{m}$ are the $m$-th row of coefficient matrices $\Theta\in\mathbb{R}^{M\times p}$ and $B\in\mathbb{R}^{M\times K}$, respectively. Without loss of generality, we assume rank$(B)=K<M$ since if $B$ is not of full column rank we can always reduce the dimensions of $Z$ such that the full column rank condition is met. Finally, we assume $Y_{m}$ and $Y_{m^{\prime}}$ are independent given $X$ and $Z$ for $m^{\prime}\neq m$.

To characterize the effect of hidden variables, we assume the following factor model (Bai, 2003; Fan et al., 2013, 2008) that relates $X$, the observed variables, to $Z$, the hidden variables:

[TABLE]

where the noise term $W\in\mathbb{R}^{p}$ has zero mean and is independent of $Z$, and $A\in\mathbb{R}^{p\times K}$ is a matrix of unknown parameters. In this paper, we focus on the independent and homogeneous noise setting where $\Sigma_{W}=\mathbb{E}(WW^{T})=\tau I_{p}$, and without loss of generality, we set $\tau=1$. The proposed method can be easily extended to the dependent noise setting, provided the smallest and largest eigenvalues of $\Sigma_{W}$ are bounded from below and above by some constants.

To understand the challenge in the identifiability of $\Theta$, we first consider a special case of (1) in which $Y_{m}$ follows the linear regression model $Y_{m}=\Theta_{m}X+B_{m}Z+E_{m}$, where $E_{m}$ is the random noise. As shown in Bing et al. (2022), the model can be rewritten as $Y_{m}=(\Theta_{m}+B_{m}L)X+\epsilon_{m}$, where $L=\mathbb{E}(ZX^{T})\{\mathbb{E}(XX^{T})\}^{-1}$ is obtained by $L_{2}(P)$ projecting $Z$ onto the linear space generated by $X$ and $\epsilon_{m}=B_{m}(Z-LX)+E_{m}$. As a result, ignoring the hidden variable $Z$ and regressing $Y_{m}$ on $X$ leads to a biased estimator of $\Theta_{m}$. In fact, it is easily seen that the coefficient matrix $\Theta+BL$ can be identified via the first two moments of $(Y,X)$. To establish the identifiability of $\Theta$, a very natural idea is to separate $\Theta$ and $BL$ in the coefficient matrix $\Theta+BL$. In the literature, a commonly used identifiability assumption for $\Theta$ is $P_{B}\Theta=0$ (Bing et al., 2022; Lee et al., 2017; Wang et al., 2017), where $P_{B}=B(B^{T}B)^{-1}B^{T}\in\mathbb{R}^{M\times M}$ is the projection matrix onto the column space of $B$. Under this assumption, the two matrices $\Theta$ and $BL$ belong to two orthogonal spaces, and therefore we can identify $\Theta$ via $\Theta=P_{B}^{\perp}(\Theta+BL)$ where $P_{B}^{\perp}=I_{M}-P_{B}$, provided $P_{B}$ is identifiable. In this case, $\Theta$ can be naturally interpreted as the association between $X$ and $Y$ that cannot be explained via the hidden variables.

Nevertheless, the above analysis suffers from the following two challenges when extended to the GLM setting in (1). First, unlike linear regression, for the GLM, the parameter obtained by regressing $Y_{m}$ on $X$ does not have a simple closed form, and therefore the relationship between this parameter and the parameter of interest $\Theta$ is unclear. Following the classical literature on misspecified models (White, 1982), a routine approach is to define the pseudo-true parameter as

[TABLE]

Since in general $b(\cdot)$ is not a quadratic function, $F_{m}^{\textsc{MLE}}$ does not have a simple closed form solution, which complicates the analysis of the identifiability of $\Theta$. To overcome this challenge, we first focus on quantifying the approximation bias of $F_{m}^{\textsc{MLE}}$ locally around the target parameter $\Theta_{m}$. Following the logic similar to the proof of Theorem 1 in Section 4, we can establish that

[TABLE]

where the first term on the right hand side corresponds to the first-order bias of $F_{m}^{\textsc{MLE}}$ and $Rem_{m}$ represents the approximation error which is of a smaller order. We note that the factor model (2) plays a pivotal role in deriving the leading bias term and quantifying the magnitude of $Rem_{m}$. In contrast, under linear regression, we have $F_{m}^{\textsc{MLE}}-\Theta_{m}=B_{m}\mathbb{E}(ZX^{T})\{\mathbb{E}(XX^{T})\}^{-1}$, which does not require the factor model (2) to hold as this result is derived purely from the $L_{2}(P)$ projection of $Z$ onto the linear space of $X$.

To identify $\Theta$, our next step is to separate $\Theta_{m}$ and the first-order bias in the decomposition of $F_{m}^{\text{MLE}}$ in (4). This brings us to the second major challenge in extending to the GLM setting, which is that the first order bias term of $F^{\text{MLE}}$ no longer lives in the column space of $B$. To see this more clearly, following the analysis used in linear regression, we stack the first-order biases in (4) over $1\leq m\leq M$ into a matrix $\mathbb{E}(DBZZ^{T})A^{T}\mathbb{E}(DXX^{T})$, where $D\in\mathbb{R}^{M\times M}$ is a diagonal matrix with the $m$th entry being $b^{\prime\prime}(\Theta_{m}X+B_{m}Z)$. Ignoring the $Rem_{m}$ term, we can write $F^{\textsc{MLE}}\approx\Theta+\mathbb{E}(DBZZ^{T})A^{T}\mathbb{E}(DXX^{T})$, where the matrix $F^{\textsc{MLE}}$ is identifiable row-by-row through (3). Unfortunately, we are not able to separate $\Theta$ and $\mathbb{E}(DBZZ^{T})A^{T}\mathbb{E}(DXX^{T})$ as in the linear case since $\mathbb{E}(DBZZ^{T})A^{T}\mathbb{E}(DXX^{T})$ is no longer in the column space of $B$ due to the presence of the matrix $D$. Consequently, under the same assumption $P_{B}\Theta=0$, the non-orthogonality of $\Theta$ and $\mathbb{E}(DBZZ^{T})A^{T}\mathbb{E}(DXX^{T})$ implies

[TABLE]

which then results in $\Theta$ not being identifiable via $P_{B}^{\perp}F^{\textsc{MLE}}$.

2.2 Our proposed approach

To address the identifiability problem, our main idea is to construct a properly weighted score function of the misspecified GLM to restore the orthogonality between $\Theta$ and the corresponding first-order bias. More precisely, for each $1\leq m\leq M$, we define the parameter $F_{m}^{*}\in\mathbb{R}^{1\times p}$ as the solution of the following estimating equation:

[TABLE]

Compared to the score function from (3), the estimating equation (5) contains a denominator $b^{\prime\prime}(F_{m}^{*}X)$, which can be viewed as an approximation of the variance of $Y_{m}$ under the GLM, i.e., $\mathrm{Var}(Y_{m}|X,Z)=b^{\prime\prime}(\Theta_{m}X+B_{m}Z)$. Thus, we can also interpret (5) as an inverse variance weighted score function. Since the GLM in (3) is misspecified, in general we have $F_{m}^{*}\neq F_{m}^{\textsc{MLE}}$. The rationale behind the weighted score approach is that the first-order bias of $F_{m}^{*}$ will have a more desirable form, which will facilitate the analysis of the identifiability of $\Theta$. Indeed, Theorem 1 shows that

[TABLE]

where $Rem_{m}^{\prime}$ presents the remainder in the expansion which is of a smaller order. It is easily seen that the first-order bias on the right hand side of (6), when stacked satisfies $P_{B}^{\perp}B\leavevmode\nobreak\ \mathbb{E}(ZZ^{T})A^{T}\{\mathbb{E}(XX^{T})\}^{-1}=0$. Under the assumption $P_{B}\Theta=0$, we can show that, ignoring the $Rem_{m}^{\prime}$ term,

[TABLE]

As a result, we can asymptotically identify $\Theta$, provided $P_{B}$ is identifiable and the $Rem_{m}^{\prime}$ term is asymptotically negligible.

To justify the identifiability of $P_{B}$, we similarly define the inverse variance weighted residual as

[TABLE]

As shown in the proof of Theorem 3, we have

[TABLE]

where $\Gamma=\mathbb{E}(ZZ^{T})A^{T}\{\mathbb{E}(XX^{T})\}^{-1}$,

[TABLE]

and $Rem^{\prime\prime}$ denotes the remainder term induced by the second-order term $Rem^{\prime}$ in (6). Recall that the first term on the right hand side of (8), $\mathbb{E}({\epsilon}{\epsilon}^{T})$, is a diagonal matrix. Consider the singular value decomposition of $B=V\Lambda U^{T}$ where $V\in\mathbb{R}^{M\times K}$ and $U\in\mathbb{R}^{K\times K}$ consist of the left and right singular vectors of $B$, respectively, and $\Lambda$ is the diagonal matrix of non-increasing singular values. From (8), under the pervasiveness assumption in the factor model literature (Bai, 2003; Fan et al., 2013, 2008), we can (asymptotically) recover $V$ by applying spectral decomposition on $\mathbb{E}(\bar{\epsilon}\bar{\epsilon}^{T})$ and obtaining the first $K$ eigenvectors. Since we can verify $P_{B}=VV^{T}$, the projection matrix $P_{B}$ is identifiable.

Compared to the analysis of the identifiability of $P_{B}$ in linear regression (Bing et al., 2022), our argument differs in the following two ways. First, the decomposition in (8) is applied to the covariance of the inverse variance weighted residual, rather than the residual itself. Again, this reweighting approach is crucial to ensure that the column space of $B$ can be identified by the spectral decomposition of a proper covariance matrix. Second, the non-linear property of $b^{\prime}(\cdot)$ requires a more careful analysis of the matrix perturbation errors in (8) to apply the Davis-Kahan Theorem. In particular, the sample version of $Rem^{\prime\prime}$ in (8) and the plug-in estimators of $F^{*}_{m}$ are correlated, leading to a slower rate of convergence. To address this technical challenge, we rely on cross-fitting and data splitting to facilitate the theory.

3 Parameter Estimation and Approximate Inference: G-HIVE

Recall that given $n$ i.i.d. observations $(Y^{(i)},X^{(i)})$ of $(Y,X)$, $i=1,...,n$, our goal is to estimate and do inference on $\Theta$. In this section, we present our estimation and inference procedure called G-hive, short for Generalized - HIdden Variable adjusted Estimation. The algorithm, inspired by the identifiability of $\Theta$ detailed in Section 2.2, is summarized in Algorithm 1.

3.1 Parameter Estimation

In this algorithm, we first randomly split the data into two folds $D_{1}$ and $D_{2}$. We use the data in $D_{1}$ to estimate the pseudo-true parameter $F_{m}^{*}$ in (5). A straightforward approach is to solve the estimating equation (5) with the expectation replaced by the sample average. However, in general, solving estimating equations may lead to multiple solutions (or the solution may not even exist), which complicates practical implementation. As an alternative, we propose to estimate $F_{m}^{*}$ by maximizing the following modified quasi-likelihood function:

[TABLE]

$Q^{(D_{1})}_{m}(F_{m})$ is a valid likelihood-type function since on the population level $\mathbb{E}(\nabla Q^{(D_{1})}_{m}(F_{m}^{*}))=0$ by (5), and $\mathbb{E}(\nabla^{2}Q^{(D_{1})}_{m}(F_{m}^{*}))$ is negative definite. This implies that, on the sample level, $Q^{(D_{1})}_{m}(F_{m})$ is locally strictly concave around $F_{m}^{*}$. However, the function $Q^{(D_{1})}_{m}(F_{m})$ may not be strictly concave for all $F_{m}$, which implies the possibility of local maximizers. Thus, following the convention in the statistics literature, it is advisable to maximize $Q^{(D_{1})}_{m}(F_{m})$ from multiple initial values in practice. Finally, we note that our function $Q^{(D_{1})}_{m}(F_{m})$ is related but different from the standard quasi-likelihood (Wedderburn, 1974) defined as $\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{b^{\prime}(F_{m}X^{(i)})}\frac{Y_{m}^{(i)}-\mu}{V(\mu)}d\mu$, where $V(\cdot)$ is the variance function. Since the pseudo-true parameter $F_{m}^{*}$ is defined under a misspecified GLM, maximizing the standard quasi-likelihood does not yield a consistent estimator of $F_{m}^{*}$.

After obtaining $\hat{F}^{(D_{1})}_{m}$ for all $1\leq m\leq M$, we can construct the estimated residuals using the data in $D_{2}$. That is, for $i\in D_{2}$, we can construct

[TABLE]

We can then estimate $\mathbb{E}(\bar{\epsilon}\bar{\epsilon}^{T})$ in (8) by

[TABLE]

where $\hat{\epsilon}^{(i)}=(\hat{\epsilon}_{1}^{(i)},...,\hat{\epsilon}_{M}^{(i)})^{T}\in\mathbb{R}^{M}$. The sample splitting procedure guarantees the desired independence between $\hat{F}^{(D_{1})}_{m}$ and the data $(Y_{m}^{(i)},X^{(i)})$ in $D_{2}$, which simplifies the technical analysis of the sample covariance matrix of $\hat{\epsilon}^{(i)}$. To fully utilize the data, we can switch the role of $D_{1}$ and $D_{2}$ to construct the estimators $\hat{F}^{(D_{2})}_{m}$ and $\hat{\Sigma}^{(D_{1})}$, and eventually define

[TABLE]

Inspired by (8), we apply spectral decomposition on $\hat{\Sigma}$ to get the first $K$ eigenvectors which are arranged as columns in $\hat{V}\in\mathbb{R}^{M\times K}$ which is then used to construct $\hat{P}_{B}^{\perp}=I-\hat{V}\hat{V}^{T}$. We refer to this as the PCA step. In view of (6), we propose to remove the first-order bias of $F^{*}$ and estimate $\Theta$ with $\hat{\Theta}=\hat{P}_{B}^{\perp}\hat{F}$, where $F^{*}:=[F_{1}^{*T},...,F_{M}^{*T}]^{T}\in\mathbb{R}^{M\times K}$ and $\hat{F}:=[\hat{F}_{1}^{T},...,\hat{F}_{M}^{T}]^{T}\in\mathbb{R}^{M\times K}$.

Remark 1**.**

Since $K$, the number of eigenvectors of $\hat{\Sigma}$ to extract, is unknown in practice, the user needs to specify its value to implement the PCA step. Similar to Ahn and Horenstein (2013); Lam and Yao (2012), we consider the following eigenvalue ratio approach. In particular, we estimate $K$ by

[TABLE]

where $\hat{\lambda}_{1}\geq\hat{\lambda}_{2}\geq...$ are the eigenvalues of $\hat{\Sigma}$ and $\bar{K}$ is a pre-specified value. Similar to Lam and Yao (2012), we set $\bar{K}=\lfloor(n\wedge M)/2\rfloor$, as the rank of $\hat{\Sigma}$ is no greater than $n\wedge M$, and it is reasonable to look at the first half of the non-zero eigenvalues. The intuition of this approach is that by (8) the sample covariance matrix $\hat{\Sigma}$ should have $K$ spiked eigenvalues and therefore the eigenvalue ratio ${\hat{\lambda}_{j}}/{\hat{\lambda}_{j+1}}$ is expected to reach the maximum at $j=K$. One desired property of this approach is that it does not require any knowledge of unknown population level quantities or additional tuning parameters. The theoretical justification of (14) follows the same argument as in Bing et al. (2022). We defer further technical results to Section C of the Appendix. In simulations, we implement this data-driven choice of $\hat{K}$ in data-driven G-hive, and it is shown to yield reasonable results.

3.2 Approximate Inference

In this subsection, we consider how to construct an inference procedure for $\Theta$. Recall that following the analysis in Section 2.2, we can only argue $\Theta$ is asymptotically identifiable. To study the inferential property for $\Theta$, we have to characterize the asymptotic bias in the identification of $\Theta$. In view of (6), Theorem 1 below shows that the asymptotic bias satisfies

[TABLE]

where $(P_{B}^{\perp}F^{*})_{m}$ is the $m$th row of $P_{B}^{\perp}F^{*}$. An important consequence of the above result is that the asymptotic bias may dominate the estimation error, making inference on $\Theta$ difficult or even infeasible.

To explain the details, we focus on the estimator $\hat{F}$. The same argument applies to $\hat{\Theta}$ as well. Theorem 2 below shows that $\hat{F}$ is asymptotically linear. That is, under some conditions, for any $1\leq m\leq M$ and $1\leq j\leq p$,

[TABLE]

where $e_{j}$ is a unit basis vector with the $j$th entry being 1 and 0 otherwise, and $G_{m}$ is defined in (18) below. Applying the central limit theorem to the first term on the right hand side and using (15), we have

[TABLE]

for some $\sigma^{2}_{mj}>0$. When $p=o(n)$, which is the regime considered in this work, the asymptotic bias dominates the stochastic error of the estimator $\hat{F}_{mj}$ so that the confidence interval based on the limiting distribution of the estimator $\hat{F}_{mj}$ does not yield the desired coverage probability for $\Theta_{mj}$. Unlike linear regression with hidden variables, the presence of asymptotic bias makes inference in GLMs substantially more challenging.

To overcome this difficulty, we propose the following approximate inference framework. Specifically, we shift the parameter of interest from $\Theta$ to $F^{*}$ or $P_{B}^{\perp}F^{*}$, where, by (15), $F^{*}$ can be viewed as the first-order approximation of $\Theta$ and $P_{B}^{\perp}F^{*}$ the second-order approximation. As a result, we refer to inference on $F^{*}$ and $P_{B}^{\perp}F^{*}$ as first-order approximate inference and second-order approximate inference, respectively. Indeed, first-order approximate inference (i.e., inference on $F^{*}$) is immediately available from the result in (16) or more generally, from Theorem 2. However, it is reasonable to expect that inference on $P_{B}^{\perp}F^{*}$ would serve as a more accurate surrogate and provide more information on $\Theta$ compared to the first-order inference method. Thus, in this work we focus on the following second-order approximate inference method for uncertainty quantification.

Assume that the parameter of interest is defined as $u^{T}(P_{B}^{\perp}F^{*})v$, where $u\in\mathbb{R}^{M}$ and $v\in\mathbb{R}^{p}$ are known vectors satisfying $||u||_{2}=||v||_{2}=1$. We thus allow for the inference on arbitrary linear combinations of parameters in our approach. Define a weighted covariance matrix as

[TABLE]

which corresponds to the expected Hessian of the modified quasi-likelihood function, where

[TABLE]

By Theorem 5, under certain conditions, we can show that

[TABLE]

where $s_{n}^{2}=\sum_{i=1}^{n}\mathbb{E}(u^{T}P_{B}^{\perp}h^{(i)})^{2}$ with $h^{(i)}=(h^{(i)}_{1},...,h^{(i)}_{M})^{T}$ and $h^{(i)}_{m}=\bar{\epsilon}^{(i)}_{m}v^{T}G^{-1}_{m}X^{(i)}$. In addition, define $\hat{h}^{(i)}=(\hat{h}^{(i)}_{1},...,\hat{h}^{(i)}_{M})^{T}$, where $\hat{h}^{(i)}_{m}=\hat{\epsilon}^{(i)}_{m}v^{T}\hat{G}^{-1}_{m}X^{(i)}$ and

[TABLE]

is an estimate of $G_{m}$. Theorem 5 further shows that the asymptotic variance $s_{n}^{2}/n$ can be consistently estimated by $\hat{s}_{n}^{2}/n$, where $\hat{s}_{n}^{2}=\sum_{i=1}^{n}(u^{T}\hat{P}_{B}^{\perp}\hat{h}^{(i)})^{2}$. As a result, the $(1-\alpha)\%$ confidence interval for $u^{T}(P_{B}^{\perp}F^{*})v$ is given by $(u^{T}\hat{\Theta}v-q_{1-\alpha/2}\hat{s}_{n}/\sqrt{n},u^{T}\hat{\Theta}v+q_{1-\alpha/2}\hat{s}_{n}/\sqrt{n})$, where $q_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of a standard normal distribution. Finally, we note that while the proposed confidence interval yields the desired coverage probability for $u^{T}(P_{B}^{\perp}F^{*})v$ rather than $u^{T}\Theta v$, we expect the proposed approximate inference framework to offer a valuable toolbox to quantify the uncertainty of estimating $\Theta$, and to provide useful information for inferring the magnitude of $\Theta$ in practice. This is confirmed in our simulation studies.

4 Statistical Guarantees

We use $\Sigma_{X}=\mathbb{E}(XX^{T})$ and $\Sigma_{Z}=\mathbb{E}(ZZ^{T})$ to denote the covariance matrices of $X$ and $Z$, respectively, and throughout the paper we consider the asymptotic setting of $p,M,K\rightarrow\infty$ as $n\rightarrow\infty$.

Assumption 1**.**

(Identifiability Assumption). Assume that $P_{B}\Theta=0$.

Assumption 2**.**

(Tail Assumption). Assume that $W$ and $Z$ are sub-Gaussian vectors with bounded sub-Gaussian norm $\sigma_{W}^{2}$ and $\sigma_{Z}^{2}$, respectively. Given $X$ and $Z$, the error $Y_{m}-b^{\prime}(\Theta_{m}X+B_{m}Z)$ is sub-exponential with bounded sub-Exponential norm $\sigma^{2}_{\epsilon,\max}$ for $1\leq m\leq M$, and \underset{1\leq i\leq n}{\max}\leavevmode\nobreak\ \underset{1\leq j\leq p}{\max}\leavevmode\nobreak\ \big{|}X^{(i)}_{j}\big{|}\leavevmode\nobreak\ \leq\leavevmode\nobreak\ C_{0} for some constant $C_{0}>0$.

Assumption 3**.**

(GLM Assumption). There exist some constants $C_{1},C_{2},C_{3}>0$, such that $C_{1}\leq b^{\prime\prime}(t)\leq C_{2}$, and $|b^{\prime}(t)|$, $|b^{\prime\prime\prime}(t)|$, $|b^{\prime\prime\prime\prime}(t)|$ are all upper bounded by $C_{3}$.

Assumption 4**.**

(Factor Model Assumption).

$\kappa_{A,1}\cdot p\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \lambda_{k}(A^{T}A)\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \kappa_{A,2}\cdot p$ for some fixed constants $\kappa_{A,1},\kappa_{A,2}>0$ and all $1\leq k\leq K$.

$\kappa_{B,1}\cdot M\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \lambda_{k}(B^{T}B)\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \kappa_{B,2}\cdot M$ for some fixed constants $\kappa_{B,1},\kappa_{B,2}>0$ and all $1\leq k\leq K$.

$\kappa_{Z,1}\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \lambda_{k}(\Sigma_{Z})\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \kappa_{Z,2}$ for some fixed constants $\kappa_{Z,1},\kappa_{Z,2}>0$ and all $1\leq k\leq K$.

$||B_{m}||_{2}\leq\leavevmode\nobreak\ C_{4}$ for all $1\leq m\leq M$ and for some fixed constant $C_{4}>0$.

Assumption 1 ensures the identifiability of $\Theta$ as explained in Section 2.2. In the related literature, there exist alternative identifiability conditions. We defer the detailed discussions to Lee et al. (2017), Bing et al. (2022), Wang et al. (2017) and Bai (2003). Assumption 2 characterizes the tail behavior of the random vectors $W$ and $Z$ and the response variable $Y$. The sub-exponential condition for $Y_{m}-b^{\prime}(\Theta_{m}X+B_{m}Z)$ holds for most GLMs such as linear, logistic and Poisson regression. To simplify the proof, we also assume the elements in $X$ are upper bounded by a fixed constant, which can be relaxed by allowing $C_{0}$ to scale with $n$ and $p$ (e.g., $C_{0}\asymp\sqrt{\log(np)}$ for sub-Gaussian $X^{(i)}_{j}$). Assumption 3 on the higher order derivatives of $b(t)$ is standard for analyzing GLMs. Finally, Assumptions 4(a) and 4(b) are known as the pervasiveness assumption in the factor model literature (Fan et al., 2013, 2008; Chang et al., 2015). A concrete example of when it is satisfied is discussed in Section A of the Appendix. Assumptions 4(c) and 4(d) are also mild conditions for factor models.

We first present a theorem that characterizes the approximation bias of $F_{m}^{*}$ defined in (5).

Theorem 1**.**

Under Assumptions 1- 4, for $1\leq m\leq M$, there exists $F_{m}^{*}$ defined in (5) such that

[TABLE]

where $\max_{1\leq m\leq M}\|Rem_{m}^{\prime}\|_{2}=O(1/p)$. In addition, the following hold:

[TABLE]

and

[TABLE]

This further implies that

[TABLE]

This theorem provides a rigorous justification of equation (6) in Section 2.2, where the remainder term $Rem_{m}^{\prime}$ in $L_{2}$ norm is of order $O(1/p)$ uniformly over $m$, and the first-order bias is of order $O(1/\sqrt{p})$ by (22). More importantly, it shows the theoretical advantage of the PCA step in our estimation procedure as it reduces the inherent bias that occurs from model misspecification due to the hidden variables. As seen in (23), using $F^{*}$ as a proxy of $\Theta$ inevitably incurs the approximation bias (in terms of the $L_{2}$ error per response) with rate $O(1/\sqrt{p})$. However, by projecting $F^{*}$ to the orthogonal space of $B$ via the PCA step, we can reduce the approximation bias to have a faster rate of $O(1/p)$.

The next theorem provides the bound for the stochastic error \big{|}\big{|}\hat{F}_{m}-F_{m}^{*}\big{|}\big{|}_{2} and the asymptotic linear approximation of $\hat{F}_{m}-F_{m}^{*}$ uniformly over $1\leq m\leq M$.

Theorem 2**.**

Under Assumptions 1- 4, $p\sqrt{\frac{\log(p\vee M)}{n}}\log(M\vee n)=o(1)$ and $\{\log(M\vee p)\}^{3}=O(n)$, there exists a local maximizer $\hat{F}^{(D_{j})}_{m}$ of $Q^{(D_{j})}_{m}(F_{m})$ for $1\leq m\leq M$ and $j\in\{1,2\}$, such that

[TABLE]

In addition, for any $v\in\mathbb{R}^{p}$ with $\|v\|_{2}=1$, we have

[TABLE]

where $G_{m}$ is defined in (18) and \max_{1\leq m\leq M}|Rem^{\prime\prime}_{m}|=O_{p}\Big{(}p^{3/2}\cdot\frac{\log(p\vee M)\log(n\vee M)}{n}\Big{)}.

Recall that the modified quasi-likelihood $Q^{(D_{j})}_{m}(F_{m})$ is non-concave and may have multiple local solutions. This theorem only applies to some local maximizer of $Q^{(D_{j})}_{m}(F_{m})$. While the rate of convergence obtained in (24) agrees with the existing literature on M-estimation with increasing dimension (Portnoy, 1984), a unique challenge that had to be overcome is the fact that the covariance matrix $\Sigma_{X}$ has spiked eigenvalues, which is implied by the hidden variable model (2) and the pervasiveness assumption in Assumption 4. This required a more delicate analysis to control the perturbation of the Hessian matrix $\nabla^{2}Q^{(D_{j})}_{m}(F_{m})$ around $F_{m}^{*}$.

Combining the approximation error in Theorem 1 and the stochastic error in Theorem 2, we obtain

[TABLE]

where in the last step we notice that the error bound is dominated by the approximation error under the condition $p\sqrt{\frac{\log(p\vee M)}{n}}\log(M\vee n)=o(1)$ in Theorem 2.

Finally, the asymptotic linear expansion of $(\hat{F}_{m}-F_{m}^{*})v$ in (25) shows that (16) holds under the condition $p^{3/2}\cdot\frac{\log(p\vee M)\log(n\vee M)}{\sqrt{n}}=o(1)$, which validates first-order approximate inference. Since inference on $F^{*}$ is not the main focus of this work, we do not pursue further results along this line.

The next theorem provides the rate for the estimation error of $\hat{P}_{B}^{\perp}$ in the PCA step.

Theorem 3**.**

Under the assumptions in Theorem 2, if we further assume $\{\log(M\vee p)\}^{5}=O(n)$ and $K<n$, then we have

[TABLE]

The estimation error of the projection matrix consists of four terms. The first two terms are the asymptotic bias corresponding to the remainder term in the expansion (8), and the last two terms stem from the stochastic error of $\hat{P}_{B}^{\perp}$. Under the condition $p\sqrt{\frac{\log(p\vee M)}{n}}\log(M\vee n)=o(1)$ in Theorem 2, the stochastic error is not necessarily dominated by the asymptotic bias in (27). So we need to keep all four terms in (27). Under additional conditions $p=o(\sqrt{M})$ and $K=o(n)$, the estimator $\hat{P}_{B}^{\perp}$ is consistent in Frobenius norm.

The proof of Theorem 3 relies on the Davis-Kahan Theorem (Lemma 6) and the bound for $\|\hat{\Sigma}-\Sigma\|_{F}$, where $\Sigma=B(\Sigma_{Z}^{-1}+A^{T}A)^{-1}B^{T}$. Based on a more refined expansion compared to (8), we can decompose the error $\hat{\Sigma}_{mm^{\prime}}-\Sigma_{mm^{\prime}}$ into pairwise interactions of 6 error terms (21 terms in total), where we further need to distinguish the analysis for the diagonal term $m=m^{\prime}$ and the off-diagonal term $m\neq m^{\prime}$. The resulting proof is much more technical than that for linear regression. In particular, we apply sample splitting to decorrelate the error terms in the expansion of $\hat{P}_{B}^{\perp}$, leading to a faster rate of convergence for some of the error terms.

The next theorem provides the rate for the estimation error of our final estimator $\hat{\Theta}$.

Theorem 4**.**

Under the same assumptions as in Theorem 3, we have

[TABLE]

where

[TABLE]

This theorem shows that the per-response $L_{2}$ estimation error of $\hat{\Theta}$ is bounded by three terms $Err_{1},Err_{2}$ and $Err_{3}$, where $Err_{1}$ is inherent from the estimation error of $\hat{P}_{B}$ in Theorem 3, $Err_{2}$ comes from the estimation error of $\hat{F}$ in Theorem 2, and $Err_{3}$ corresponds to the approximation error of $P_{B}^{\perp}F^{*}$ in Theorem 1. When $p=o(\sqrt{M}),\leavevmode\nobreak\ K=o(n)$, $p\sqrt{\frac{\log(p\vee M)}{n}}\log(M\vee n)=o(1)$, and $||F^{*}||_{\text{op}}/\sqrt{M}=O(1)$, as $p,M,n,K\rightarrow\infty$, the estimator $\hat{\Theta}$ is consistent in the sense that $\|\hat{\Theta}-\Theta\|_{F}/\sqrt{M}=o(1)$.

Remark 2**.**

To have a more refined comparison with the estimation error of $\hat{F}$ in (26), we further assume $||\Theta||_{\text{op}}=O(1)$. Under this assumption, since by Theorem 1 we have ||F^{*}||_{\text{op}}\leq||\Theta||_{\text{op}}+\|F^{*}-\Theta\|_{F}=O\big{(}1+\sqrt{M/p}\big{)}, the rate of the estimator $\hat{\Theta}$ in Theorem 4 reduces to

[TABLE]

Assuming $M\asymp p^{\alpha}$ and $K\asymp n^{\beta}$ for some positive constants $\alpha\geq 1$ and $\beta\leq 1$, the rate can be simplified to

[TABLE]

Provided $\alpha>2$ and $\beta<1$, the rate of our estimator $\hat{\Theta}$ in (29) is faster than the rate of $\hat{F}$ in (26), which justifies the theoretical benefit of our proposed method over the naive MLE approach that ignores hidden variables.

Recall the notation in Section 3.2: $G_{m}$ and $\hat{G}_{m}$ are defined in (18) and (20), $s_{n}^{2}=\sum_{i=1}^{n}\mathbb{E}(u^{T}P_{B}^{\perp}h^{(i)})^{2}$ with $h^{(i)}=(h^{(i)}_{1},...,h^{(i)}_{M})^{T}$ and $h^{(i)}_{m}=\bar{\epsilon}^{(i)}_{m}v^{T}G^{-1}_{m}X^{(i)}$, and $\hat{s}_{n}^{2}=\sum_{i=1}^{n}(u^{T}\hat{P}_{B}^{\perp}\hat{h}^{(i)})^{2}$ with $\hat{h}^{(i)}=(\hat{h}^{(i)}_{1},...,\hat{h}^{(i)}_{M})^{T}$ and $\hat{h}^{(i)}_{m}=\hat{\epsilon}^{(i)}_{m}v^{T}\hat{G}^{-1}_{m}X^{(i)}$. Finally, we establish the limiting distribution of the estimator $u^{T}\hat{\Theta}v$ in the following theorem.

Theorem 5**.**

Under the same assumptions as in Theorem 3, if we further assume $K=o(n)$, $p=o(\sqrt{M})$ and $\mathbb{E}(u^{T}P_{B}^{\perp}h^{(i)})^{2}\geq C$ for some constant $C>0$, then for any $u\in\mathbb{R}^{M}$ and $v\in\mathbb{R}^{p}$ with $\|u\|_{2}=\|v\|_{2}=1$,

[TABLE]

for some constant $C^{\prime}>0$, where $\Phi(\cdot)$ is the c.d.f of a standard normal distribution, and the three error terms in (30) are given by

[TABLE]

[TABLE]

and

[TABLE]

with $\|u\|_{1}\leq R_{n}$ for some $R_{n}>0$. Finally, the asymptotic variance $s_{n}^{2}/n$ can be consistently estimated by $\hat{s}_{n}^{2}/n$, i.e.,

[TABLE]

We characterize the accuracy of the Gaussian approximation of $u^{T}\hat{\Theta}v$ via the Berry–Esseen bound (30), where $\delta_{1}$ corresponds to the Gaussian approximation error of $\hat{F}$ which comes from the asymptotic linear expansion (25) in Theorem 2, $\delta_{2}$ corresponds to the estimation error of $\hat{P}_{B}$, and $\delta_{3}$ corresponds to the product error of $\hat{F}$ and $\hat{P}_{B}$. To simplify the Berry–Esseen bound (30), we further assume that $R_{n}=O(1)$, $K=O(1)$ and $\|F^{*}v\|_{2}=O(1)$. Then the Berry–Esseen bound (30) reduces to

[TABLE]

Assuming $M\asymp n^{r_{M}}$ and $p\asymp n^{r_{p}}$ for some positive constants $r_{M}$ and $r_{p}$, the above bound goes to 0 when $r_{p}<1/3$ and $r_{M}>(1-r_{p})\vee(\frac{1}{2}+r_{p})$ hold. The condition $r_{p}<1/3$ on the number of covariates $p$ is comparable to the requirement on the dimensionality for M-estimation or the intrinsic dimensionality for sparse GLMs. The condition $r_{M}>(1-r_{p})\vee(\frac{1}{2}+r_{p})$ is unique in our hidden variable model, which implies that the number of responses $M$ needs to be large enough such that we can borrow information from different responses to better estimate $P_{B}$. It can also be viewed as a kind of blessing of dimensionality, as the error (32) generally decreases with larger $M$.

Finally, as long as the asymptotic variance $s_{n}^{2}/n$ can be consistently estimated by $\hat{s}_{n}^{2}/n$ shown by (31), the second-order approximate confidence interval proposed in Section 3.2 yields the desired coverage probability for $u^{T}(P_{B}^{\perp}F^{*})v$.

5 Simulation Results

Here we present our simulation results, which can be divided into three categories: the approximation bias of $F^{*}$ and $P_{B}^{\perp}F^{*}$, the estimation error of $\hat{\Theta}$, and statistical inference of $\hat{\Theta}$. We first present the data generating mechanism, and then discuss each result in the above categories.

5.1 The Data Generating Mechanism

To satisfy Assumption 4, we first generate $pK$ i.i.d. $N(0,1)$ random variables to construct the matrix $A\in\mathbb{R}^{p\times K}$ and normalize each of the $p$ rows to have an $L_{2}$ norm of 1. The usage of Gaussian random variables to satisfy Assumption 4 is formally justified in Section A in the Appendix. We then generate $n$ i.i.d random vectors $Z\sim N(\mathbf{0}_{K},\Sigma_{Z})$, where $\Sigma_{Z}$ is a circulant matrix with $1$’s on the diagonal and a decay rate of $-0.5$. Each component of $W$ is drawn independently from $N(0,1)$ and we compute $X=AZ+W$ to generate the observed covariate matrix $X\in\mathbb{R}^{p\times n}$. Similarly, each component of the hidden coefficient matrix $B\in\mathbb{R}^{M\times K}$ is i.i.d. $N(0,1)$ and we normalize each of the $M$ rows to have an $L_{2}$ norm of 1. Then we multiply $B$ by a positive scalar $\eta\in\mathbb{R}$, where $\eta$ is a parameter we can tweak to alter the influence of the hidden variables. A larger $\eta$ value corresponds to larger confounding from the hidden variables $Z$. Each element of the coefficient matrix $\Theta\in\mathbb{R}^{M\times p}$ is i.i.d. $N(0,1)$ and we normalize each row to have an $L_{2}$ norm of 1. We project $\Theta$ onto the orthogonal column space of $B$ to satisfy Assumption 1 and get the final value for $\Theta$. Lastly, we generate $n$ i.i.d. copies of $Y_{m}\in\{0,1\}$ from a Bernoulli distribution with $\mathbb{P}(Y_{m}=1|X,Z)=\exp(\Theta_{m}X+B_{m}Z)/[1+\exp(\Theta_{m}X+B_{m}Z)]$ for each $1\leq m\leq M$. Throughout the simulations, we set $K=3$.

5.2 The Methods

We consider three variants of our G-hive algorithm. data driven G-hive is the proposed method that uses a data driven estimate of $K$ mentioned in Remark 1. The other two are oracle type estimators, Oracle(K) G-hive, corresponding to the algorithm with the true value of $K$ given and Oracle(P) G-hive, the algorithm with the true projection matrix $P_{B}$ given. These two oracle type estimators illustrate the impact of estimating $K$ with $\hat{K}$ and the impact of estimating $P_{B}^{\perp}$ with $\hat{P}_{B}^{\perp}$ in our method. The baseline method we compared against was the naive maximum likelihood estimator that ignores the hidden variables. This is denoted as naive mle and was implemented with the glm function in R.

5.3 Evaluating the Approximation Bias

As seen from Theorem 1, the first-order approximation bias \big{|}\big{|}F^{*}-\Theta\big{|}\big{|}_{F}/\sqrt{M} and the projected second-order approximation bias \big{|}\big{|}P_{B}^{\perp}F^{*}-\Theta\big{|}\big{|}_{F}/\sqrt{M} decays with $p$ with order $O(1/\sqrt{p})$ and $O(1/p)$, respectively. This characterizes the inherent and unavoidable gap between the \saytrue parameter in the misspecified GLM and the true parameter in the correctly specified GLM. We set $M=3$, $\eta=10$, and the results were averaged over $r=20$ repetitions. As obtaining the \saytrue $F^{*}$ amounts to finding the solution to the estimating equation (5) for each $1\leq m\leq M$ which does not have a closed form solution, we instead compute the solution to the modified quasi-likelihood function in (10) with an extremely large $n=2\times 10^{5}$ via the Monte Carlo approach. The results are shown in the left graph in Figure 1. It is apparent that as $p$ increases, both forms of the approximation bias indeed decay with the projected bias being much smaller than the non-projected bias. Surprisingly, the projected approximation bias is very close to 0, even if $p$ is as small as 3. This confirms the theory in Theorem 1 and also validates the projection step in our G-hive method.

5.4 Evaluating the Estimation Error of $\hat{\Theta}$

Next we evaluate the estimation error of $\hat{\Theta}$ in three scenarios: (1) when we vary the level of influence of the hidden variables through the magnitude of $\eta$, (2) when we increase the sample size $n$, and (3) when we increase the dimension of the response $M$. For (1) we vary $\eta$ to take values in $\{1,2,...,8\}$ and have the setting of $p=M=4,\leavevmode\nobreak\ K=3,\leavevmode\nobreak\ n=100$ averaged over $r=500$ repetitions. The results are shown in the right graph in Figure 1. Recall that the larger the $\eta$ value, the larger the effect of the hidden variables and thus, a more challenging simulation setting. It is not surprising that the naive mle method performs gradually worse, as the magnitude of $\eta$ increases, the MLE is obtained with respect to a model that is becoming more and more misspecified. In contrast, the three G-hive based methods are more robust to the effect of $\eta$, which verifies that the signature projection step in G-hive mitigates the effects of misalignment between the misspecified model and the true model. Lastly, it is reasonable to see an increase in performance as we supply more model information to the estimators such as the true $K$ value or the true projection matrix $P_{B}$.

For (2), to verify the consistency of our estimator, we gradually increase $n$ to take values from $100$ to $400$ in increments of $50$ and average over $r=200$ repetitions with the model setting being identical otherwise ($M=p=4,\leavevmode\nobreak\ K=3,\leavevmode\nobreak\ \eta=4$). The results are shown in the left graph of Figure 2, and it is apparent that all four methods show an improvement in performance as $n$ increases.

For (3), to explore performance in high-dimensional response settings, we vary $M$ to take values in $\{4,8,12,16,20\}$ and have the setting of $p=4,\leavevmode\nobreak\ K=3,\leavevmode\nobreak\ n=200$ averaged over $r=100$ repetitions. The results are shown in the right graph of Figure 2 and they indicate that as $M$ grows, the estimation error slightly increases for all four methods. This coincides with the estimation error of our estimator discussed in Remark 2 in that the term $\sqrt{p\log(p\vee M)/n}$ in (28) scales with $\sqrt{\log M}$. The three G-hive estimators uniformly outperforming the naive mle method shows that G-hive is viable and competitive for a wide range of response dimension values.

5.5 Evaluating the asymptotic normality of $\hat{\Theta}$

Lastly we evaluate the asymptotic normality of $\hat{\Theta}$ with our proposed data driven G-hive. We are in the similar setting of $M=p=4,\leavevmode\nobreak\ K=3,\leavevmode\nobreak\ \eta=4$ and have results for $n=40$ and $n=70$ averaged over $r=100$ repetitions with $\alpha=0.05$. While Theorem 5 only provides coverage guarantees for $P^{\perp}_{B}F^{*}$, we also include the coverage probability pertaining to the true $\Theta$. Additionally, while Theorem 5 provides inference guarantees for $u^{T}P^{\perp}_{B}F^{*}v$ for any real vectors $u\in\mathbb{R}^{M},\leavevmode\nobreak\ v\in\mathbb{R}^{p}$, for simplicity, we fix $u=v=[1,0,0,0]^{T}$ corresponding to $(P^{\perp}_{B}{F^{*}})_{11}$. We include in Table 1 the estimated standard error $\hat{s}_{n}/\sqrt{n}$, the confidence interval lengths, and the coverage probabilities for both our data driven G-hive method and the naive mle method implemented with the glm function in R. For naive mle, the constructed confidence intervals are overly narrow which is to be expected for Wald confidence intervals for misspecified models. Coverage probabilities worsen as the sample size increases, which is to be expected as model misspecification bias is left unaddressed. Thus, ignoring the hidden variables when fitting the GLM yields misleading inference results. On the other hand, for data driven G-hive, the coverage probabilities closely match the $1-\alpha$ level, showing the validity of the constructed confidence interval. It is interesting to see that the coverage probability for the true $\Theta$ is close to $1-\alpha$ as well. This implies that while the approximation bias dominates the stochastic error of the estimator and prevents us from having standard inference results for $\Theta$ as shown in (17), in practice, for large enough $n$, the inference results for $P^{\perp}_{B}F^{*}$ can be used as a proxy for $\Theta$.

All in all, these simulation results highlight the soundness of the theory and demonstrate the feasibility and advantages of our G-hive method.

6 Real Data Analysis

We apply our data driven G-hive procedure to a dataset from Chicco and Rovelli (2019) regarding mesothelioma, a type of lung cancer. Specifically, this dataset consists of real electronic health records on 324 patients in Turkey of which 96 are diagnosed with mesothelioma and 228 are not. The dataset includes 33 explanatory variables including age, platelet count, white blood cell count, etc., but in order to facilitate the analysis, we included the explanatory variables that were shown to be meaningful in terms of reducing the mean square error (MSE) in Chicco and Rovelli (2019). We removed categorical variables and variables that were highly correlated (a correlation coefficient $\geq 0.91$). We further removed 2 explanatory variables, \sayasbestos exposure and \sayduration of asbestos exposure and considered them to be hidden confounders and treated these variables as our unobserved $Z$ variables. In the end, we ran the data analysis on 9 continuous covariates.

It is known in the medical and biology literature that asbestos, a term applied to mineral species that occur in fibrous forms, can cause chronic inflammation (Committee on Asbestos: Selected Health Effects, 2006). It is also well known that an inflammatory response causes changes in the components of our blood such as white blood cell count, etc. Thus, it is straightforward to see that asbestos related explanatory variables closely affect other explanatory variables in the dataset (white blood cell count (WCC), etc.). Also, according to Chicco and Rovelli (2019), long exposure to asbestos makes mesothelioma very likely. Thus, asbestos related variables affect the response variable \saydiagnosis of mesothelioma as well. Hence, if we remove asbestos related variables from the dataset, we can consider them hidden confounders that affect both the observed covariates and the response variables and as a result align the real dataset with the model setup we have for our method, G-hive.

To construct a multivariate response data structure, we include three symptom related response variables from the dataset (\saychest ache, \saydyspnoea, \saypatient’s ability to perform normal tasks) along with the main response of interest, \saydiagnosis of mesothelioma. Thus, we have $M=4,\leavevmode\nobreak\ p=9,\leavevmode\nobreak\ n=324$ for our setting. All of the explanatory variables were standardized to have 0 mean and unit variance prior to the data analysis. The resulting $\widehat{\Theta}_{1\cdot}$, i.e. the coefficients pertaining to the diagnosis of mesothelioma are shown below:

As it is impossible to know the ground truth, we rely on the results provided in Chicco and Rovelli (2019) to gauge the accuracy of our method. The authors in Chicco and Rovelli (2019) conclude \saylung side and \sayplatelet count to be the two most important variables in classifying whether a patient has mesothelioma or not, and this is consistent with the findings with our method, even in the setting of having the asbestos related variables considered hidden confounders and removed. In terms of magnitude, the coefficient values pertaining to \saylung side and \sayplatelet count are larger than the other features. Since all of the features were standardized beforehand, this is a good indication that our method produces reasonable results. Additionally, in Chicco and Rovelli (2019), the authors claim that there is a positive correlation between \saylung side and the \saydiagnosis of mesothelioma, while there is a negative correlation between \sayplatelet count and the \saydiagnosis of mesothelioma. Our results are aligned with this fact from the literature as well, since the corresponding estimates given in Table 2 are positive and negative, respectively. While the naive mle method showed similar results, the biggest difference was in the magnitude of the most important covariate, \sayPlatelet, for which our G-hive method better represented the stark effect. Thus, our results appear to be in line with the current medical literature, even in the presence of hidden confounders.

We also applied our G-hive procedure to analyze another NHANES dataset (Centers for Disease Control and Prevention and National Center for Health Statistics (2018)). We focused on the general estimation ability of G-hive and also highlighted the effect of hidden variables in the context of confounding. The detailed results and discussion are deferred to Appendix E.

7 Discussion

In this paper we introduced G-hive, a unified framework and implementable pipeline for estimation and approximate inference in multivariate response GLMs with hidden variables that combines a novel bias correcting step with reweighted estimating equations and a spectral decomposition based projection step. More specifically, we define a novel pseudo-parameter $F^{*}$ via an inverse variance reweighted score, then remove its leading bias term by projecting onto the orthogonal complement of the latent factor column space via PCA on the covariance matrix of reweighted residuals. Theoretically, we derive the convergence rates of the first and second-order approximations, $F^{*}$ and $P_{B}^{\perp}F^{*}$, to the true parameter, $\Theta$. We also establish convergence rates for the estimation error of $\hat{F}$, $\hat{P}_{B}$, and the final proposed estimator, $\hat{\Theta}=\hat{P}_{B}^{\perp}\hat{F}$. A Berry–Esseen bound that leads to valid Gaussian inference for linear combinations of $P_{B}^{\perp}F^{*}$ is also derived. Empirically, our simulation results show that G-hive is much more robust to confounding compared to the baseline method in multiple $p,M,n$ settings, and these robustness and deconfounding benefits of G-hive were shown to extend to real-data analyses with lung cancer data and the NHANES dataset.

Our approach relies on standard but substantive assumptions. Directions for extending the current work include handling the general $\Sigma_{W}$ setting and the high-dimensional $p>n$ setting via regularized estimators, but we leave these topics for future study.

Appendix A Remark on Assumption 4

In this section we provide an example of when Assumption 4 is satisfied. Suppose each element in the $M\times K$ matrix $B$ is i.i.d. from a centered Gaussian distribution with a bounded variance. Without loss of generality, let’s assume the variance $\sigma^{2}=1$. By standard high probability bounds for Gaussian covariance estimation (Vershynin (2018)) and Weyl’s inequality (Horn and Johnson (1994)), it can be shown that the smallest eigenvalue of $B^{T}B$ is on the order of $M$. The same logic can be applied to $A$ to show that its smallest eigenvalue is on the order of $p$ with high probability. Standard Gaussian covariance estimation bounds give us with probability $\geq 1-2e^{-t}$ that

[TABLE]

for some constant $C>0$ that does not depend on $M$. This holds because we can regard the $M$ columns in $B^{T}$ as independent realizations of $N(0,I_{K})$. Then we have from Weyl’s inequality that

[TABLE]

where the last line holds with probability $\geq 1-2e^{-t}$. Choosing $t=\log M$, it is apparent that

[TABLE]

Appendix B Collection of Proofs

Lemma 1**.**

Under Assumptions 1- 4, we have

[TABLE]

for some fixed constant $c>0$, where $\tilde{Z}:=\Sigma_{Z}A^{T}(A\Sigma_{Z}A^{T}+I_{p})^{-1}X$. In addition, $\Sigma_{X}^{-1/2}X$ is a centered sub-Gaussian random vector with bounded sub-Gaussian norm.

Proof.

From the definition of $\tilde{Z}:=\Sigma_{Z}A^{T}(A\Sigma_{Z}A^{T}+I_{p})^{-1}X$ and the latent factor model $X=AZ+W$, we have the following:

[TABLE]

where the last line follows from the block matrix inverse formula. Since $W$ and $Z$ are sub-Gaussian random vectors, we merely need to compute the Euclidean norms of $B_{m}(\Sigma_{Z}^{-1}+A^{T}A)^{-1}A^{T}$ and $B_{m}(\Sigma_{Z}^{-1}+A^{T}A)^{-1}\Sigma_{Z}^{-1}$,

[TABLE]

Note that $\Sigma_{X}^{-1/2}X=\Sigma_{X}^{-1/2}AZ+\Sigma_{X}^{-1/2}W$. Following a similar derivation, we can show that $\|\Sigma_{X}^{-1/2}A\|_{\textrm{op}}$ and $\|\Sigma_{X}^{-1/2}\|_{\textrm{op}}$ are upper bounded by a constant. Thus, $\Sigma_{X}^{-1/2}X$ is a centered sub-Gaussian vector. This completes the proof. ∎

B.1 Proof for Theorem 1

Proof.

For ease of notation, for this proof, we denote $F_{m}^{*}$ as $F_{m}$. Recall from (9) that

[TABLE]

Thus, $Y_{m}=b^{\prime}(\Theta_{m}X+B_{m}Z)+\epsilon_{m}\cdot b^{\prime\prime}(\Theta_{m}X+B_{m}Z)$. Plugging this into (5) which is

[TABLE]

we get

[TABLE]

The second term is [math] as $\mathbb{E}(\epsilon_{m}|X,Z)=0$. By Taylor expansion we have

[TABLE]

for some intermediate $\delta_{m}$ between $F_{m}X$ and $\Theta_{m}X+B_{m}Z$. Rearranging we get

[TABLE]

[TABLE]

We define the following set

[TABLE]

where $\gamma>0$ is to be specified later on. In the following, we will show that $F_{m}$ implicitly given by (33) belongs to $\Omega_{m}$. To this end, we plug (33) into the condition in (34) to obtain the $\gamma$ that will satisfy this inequality.

[TABLE]

where $\tilde{Z}:=\Sigma_{Z}A^{T}(A\Sigma_{Z}A^{T}+I_{p})^{-1}X$. Using the inequality $(a+b)^{2}\leq 2a^{2}+2b^{2}$ twice, we get:

[TABLE]

and

[TABLE]

The first term in (35) can be bounded by noting the sub-Gaussian property $||B_{m}(Z-\tilde{Z})||_{\psi_{2}}\leq c/\sqrt{p}$ in Lemma 1,

[TABLE]

For the second term in (35), Lemma 1 implies $\Sigma_{X}^{-1/2}X$ is a centered, sub-Gaussian random vector, where $\Sigma_{X}^{-1/2}X$ has a bounded $\psi_{2}$-norm, i.e., that $||v^{T}(\Sigma_{X}^{-1/2}X)||_{\psi_{2}}\leq c_{x}$ for any $||v||_{2}=1$ for some fixed constant $c_{x}>0$. Then, following the definition of the sub-Gaussian norm, like in (36), we derive:

[TABLE]

For ease of notation, let us denote

[TABLE]

Then we have the following algebraic derivation:

[TABLE]

where the last step follows from (36). The first term in (38) can be bounded by:

[TABLE]

where the third line follows from the Cauchy-Schwartz inequality and the last step follows from (37).

Thus, combining these results with (35), we get the following

[TABLE]

By setting $\gamma^{4}=C/p^{2}$ for some constant $C$ sufficiently large (where $C$ does not depend on $m$), we have $\frac{128c^{4}}{p^{2}}+(\frac{c_{x}^{2}\cdot C_{3}}{C_{1}})^{4}\cdot\gamma^{8}\leq\gamma^{4}$, which implies $F_{m}^{*}\in\Omega_{m}$.

In the following, we will bound $||F_{m}-\Theta_{m}||_{2}^{2}$. Note from line (33), we have

[TABLE]

Note that from Assumption 4, the first term can be upper bounded by

[TABLE]

The second term can be upper bounded by

[TABLE]

where in the last step Assumption 3 implies \Big{|}\frac{b^{\prime\prime\prime}(\delta_{m})}{2b^{\prime\prime}(F_{m}X)}\Big{|} is bounded, \mathbb{E}[\big{(}\Theta_{m}X+B_{m}Z-F_{m}X\big{)}^{4}]\lesssim 1/p^{2} due to $F_{m}^{*}\in\Omega_{m}$, and we also apply Assumption 4. Thus, we have $||F_{m}-\Theta_{m}||_{2}^{2}\lesssim(1/p)$ uniformly over $m$, and $Rem^{\prime}$ (the term in (39)) satisfies $\max_{1\leq m\leq M}\|Rem_{m}^{\prime}\|^{2}_{2}=O(1/p^{2})$.

Compiling this into a matrix with $M$ rows, we get the following:

[TABLE]

It is trivial to note that

[TABLE]

since we get the rates for $||F_{m}-\Theta_{m}||_{2}^{2}$ uniformly over $1\leq m\leq M$. To get the last result in the theorem, we multiply both sides of (40) by $P_{B}^{\perp}$. This eliminates the first term as the $P_{B}^{\perp}B=0$. Thus, we are left with

[TABLE]

where the last line follows from the derivation in (39). This concludes the proof. ∎

B.2 Proof of Theorem 2

Since sample splitting does not change the proof, for simplicity we omit sampling splitting in the proof and define the estimator $\hat{F}_{m}$ as the local maximizer of $Q_{m}(F)$ using all $n$ data points. Let $L_{m}(F)=-Q_{m}(F)$. Our goal is to show that there exists a local minimizer $\hat{\Delta}_{m}$ of $L_{m}(F_{m}^{*}+\Delta)$ such that $\hat{\Delta}_{m}\in\mathcal{C}$ for all $1\leq m\leq M$, where $\mathcal{C}=\{\Delta\in\mathbb{R}^{p}:\|\Delta\|_{2}\leq r\}$ and $r=C\sqrt{p\log(M\vee p)/n}$ for some constant $C$ large enough. To this end, it suffices to show that the event

[TABLE]

holds with probability tending to 1, where $\partial\mathcal{C}=\{\Delta\in\mathbb{R}^{p}:\|\Delta\|_{2}=r\}$. Applying the mean value theorem, we have for any $\Delta\in\partial\mathcal{C}$,

[TABLE]

under the following two events,

[TABLE]

where $C^{\prime}$ and $C^{\prime\prime}$ are constants. Provided $C>2C^{\prime}/C^{\prime\prime}$, we obtain that $L_{m}(F_{m}^{*}+\Delta)-L_{m}(F_{m}^{*})>0$. In the following, we will show that $\mathbb{P}(E_{1})\rightarrow 1$ and $\mathbb{P}(E_{2})\rightarrow 1$. Recall that

[TABLE]

By definition we notice that $\nabla L_{m}(F_{m}^{*})$ has mean zero. Since $b^{\prime\prime}(F_{m}^{*}X^{(i)})\geq C_{1}$, all entries of $X^{(i)}$ are bounded, and $\|Y_{m}-b^{\prime}(F_{m}^{*}X)\|_{\psi_{1}}\leq\|Y_{m}-b^{\prime}(\Theta_{m}X+B_{m}Z)\|_{\psi_{1}}+\|b^{\prime}(\Theta_{m}X+B_{m}Z)-b^{\prime}(F_{m}^{*}X)\|_{\psi_{1}}$ is bounded, Bernstein inequality implies that \mathbb{P}\Big{(}|(\nabla L_{m}(F_{m}^{*}))_{j}|\geq C^{\prime}\sqrt{\frac{\log(M\vee p)}{n}}\Big{)}\leq(p\vee M)^{-3} for some constant $C^{\prime}$. Applying the union bound, we can show that

[TABLE]

with high probability, which further implies $\mathbb{P}(E_{1})\rightarrow 1$.

For the event $E_{2}$, let us denote $\tilde{F}_{m}=F_{m}^{*}+t\Delta$ and $\zeta^{(i)}_{m}(\tilde{F}_{m})=\frac{(Y^{(i)}_{m}-b^{\prime}(\tilde{F}_{m}X^{(i)}))b^{\prime\prime\prime}(\tilde{F}_{m}X^{(i)})}{\{b^{\prime\prime}(\tilde{F}_{m}X^{(i)})\}^{2}}$. We first normalize $\nabla^{2}L_{m}(\tilde{F}_{m})$ since $X$ has spiked eigenvalues under our factor model assumption. Specifically,

[TABLE]

since the smallest eigenvalue of $\Sigma_{X}=A\Sigma_{Z}A^{T}+I$ is no smaller than 1. For notational simplicity, we set $\nabla^{2}\tilde{L}_{m}(\tilde{F}_{m})=\Sigma_{X}^{-1/2}\nabla^{2}L_{m}(\tilde{F}_{m})\Sigma_{X}^{-1/2}$. Then we have

[TABLE]

where $\tilde{X}^{(i)}=\Sigma_{X}^{-1/2}X^{(i)}$, and

[TABLE]

By Lemma 1, $\tilde{X}^{(i)}$ is a sub-Gaussian vector, so $\tilde{X}_{j}^{(i)}\tilde{X}_{k}^{(i)}$ is sub-exponential. Therefore, due to the boundedness of $|b^{\prime\prime}|,\leavevmode\nobreak\ |b^{\prime\prime\prime}|$ and the fact that $Y_{m}^{(i)}-b^{\prime}(F_{m}^{*}X^{(i)})$ is sub-exponential from Assumption 2 and 3, $(1+\zeta^{(i)}_{m}(F^{*}_{m}))\tilde{X}_{j}^{(i)}\tilde{X}_{k}^{(i)}$ is $1/2$-sub-exponential. Lemma 10 implies that,

[TABLE]

provided $\{\log(M\vee p)\}^{3}=O(n)$, which implies

[TABLE]

with high probability. Furthermore, writing out $\zeta_{m}^{(i)}$ and $\zeta_{m}^{(i)^{\prime}}$, under Assumption 3, it is apparent that \Big{|}\zeta_{m}^{(i)^{\prime}}\Big{|} is bounded by constants and a random term, $F_{m}X^{(i)}$. Thus, $\zeta^{(i)}_{m}(F_{m})$ is Lipschitz in terms of $F_{m}X^{(i)}$. Using this fact, the Mean Value Theorem, and the rate for the maximum of a sub-exponential random variable, we get

[TABLE]

This then implies that

[TABLE]

where the second lines follows from plugging in our choice of rate for $||\Delta||_{2}=r=\sqrt{p\log(p\vee M)/n}$ and the moment bound $\mathbb{E}\|\frac{1}{n}\sum_{i=1}^{n}\tilde{X}^{(i)}\tilde{X}^{(i)T}-I_{p}\|_{\text{op}}\lesssim\|I_{p}\|_{\text{op}}\sqrt{p/n}$. For details of the latter, see the proof of Lemma 4. Combining these results, Weyl’s inequality implies

[TABLE]

uniformly over $m$. Therefore,

[TABLE]

Finally, we focus on

[TABLE]

where $\xi$ is an intermediate value between $\Theta_{m}X^{(i)}+B_{m}Z^{(i)}$ and $F^{*}_{m}X^{(i)}$. Notice that

[TABLE]

for some constant $C$. To lower bound $\lambda_{\min}(\mathbb{E}(1+\zeta^{(i)}_{m}(F^{*}_{m}))\tilde{X}^{(i)}\tilde{X}^{(i)T})$, we apply the following truncation technique. By Weyl’s inequality again, we have

[TABLE]

where $\Psi_{m}=\Theta_{m}X^{(i)}+B_{m}Z^{(i)}-F^{*}_{m}X^{(i)}$. Thus,

[TABLE]

For the second term on the right hand side of (43), we have

[TABLE]

where the second line follows from the Cauchy-Schwarz inequality, the third line is from the sub-Gaussian property of $\tilde{X}_{j}^{(i)}$ and the Cauchy-Schwarz inequality again, and the last step is due to $\|\Psi_{m}\|_{\psi_{2}}\lesssim p^{-1/2}$ implied by the proof of Theorem 1. Plugging these results into (43), we have

[TABLE]

which implies $\mathbb{P}(E_{2})\rightarrow 1$ in view of (41). As a by-product, we have

[TABLE]

The rest of the proof is to show (25). Recall that $\nabla L_{m}(\hat{F}_{m})=0$. This combined with the mean value theorem gives us that

[TABLE]

where $\tilde{\zeta}^{(i)}_{m}=\zeta^{(i)}_{m}(\tilde{F}_{m})$. Then for any $v\in\mathbb{R}^{p}$,

[TABLE]

where $\zeta^{(i)}_{m}=\zeta^{(i)}_{m}(F^{*}_{m})$. So it remains to control the two terms in (46) and (47) respectively. From the analysis of the term $I_{1}$ and event $E_{1}$ defined previously, using the identity $A^{-1}-B^{-1}=-A^{-1}(A-B)B^{-1}$, the term (46) is upper bounded by

[TABLE]

where we use (45), combined with the concentration bound for the term $I_{1}$ above and Weyl’s inequality to show that the smallest eigenvalue of $\frac{1}{n}\sum_{i=1}^{n}(1+\zeta^{(i)}_{m})X^{(i)}X^{(i)T}$ is also lower bounded by a constant. Following a similar argument, we can rewrite (47) as

[TABLE]

where $\tilde{H}=\frac{1}{n}\sum_{i=1}^{n}(1+\tilde{\zeta}^{(i)}_{m})\tilde{X}^{(i)}\tilde{X}^{(i)T}$ and $\hat{H}=\frac{1}{n}\sum_{i=1}^{n}(1+\zeta^{(i)}_{m})\tilde{X}^{(i)}\tilde{X}^{(i)T}$. Recall that from the analysis of the term $I_{2}$ above we have

[TABLE]

and from the analysis of the term $I_{1}$ , we have

[TABLE]

where the minimum eigenvalue of $\mathbb{E}\hat{H}$ is lower bounded by a constant as shown in (44). Thus, (47) is upper bounded by O_{p}\big{(}p^{3/2}\frac{\log(p\vee M)\log(n\vee M)}{n}\big{)}. This completes the proof of (25).

B.3 Proof for Theorem 3

Proof.

Recall that the three residuals are

[TABLE]

We use the upper-script $(i)$ to indicate the r.v from the $i$th sample. For simplicity, we use $\mathbb{E}_{n}\epsilon$ to denote $\frac{1}{n}\sum_{i=1}^{n}\epsilon^{(i)}$ and $\mathbb{E}_{n}(\epsilon|X)$ to denote $\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}(\epsilon^{(i)}|X^{(i)})$. Finally, due to sample splitting, for simplicity we can just equivalently assume that $\hat{F}_{m}$ is independent of $Y^{(i)}$ and $X^{(i)}$.

We also have $\hat{P}^{\perp}_{B}=I-\hat{V}\hat{V}^{T}$, $P^{\perp}_{B}=I-VV^{T}$, where $\hat{V}$ is the first $K$ eigenvectors of $\hat{\Sigma}=\frac{1}{n}\hat{\epsilon}^{(i)}\hat{\epsilon}^{(i)T}$ and $V$ is the first $K$ left singular vectors of $B$ and is also the first $K$ eigenvectors of $\mathbb{E}[B(Z-\tilde{Z})(Z-\tilde{Z})^{T}B^{T}]=B(\Sigma_{Z}^{-1}+A^{T}A)^{-1}B^{T}$. We thus recall, define, and denote $\Sigma=B(\Sigma_{Z}^{-1}+A^{T}A)^{-1}B^{T}$ and $\hat{\Sigma}=\mathbb{E}_{n}\hat{\epsilon}^{\otimes 2}$ accordingly.

By the Davis-Kahan theorem from Lemma 6, we get the following inequality:

[TABLE]

where $O$ is some orthogonal matrix. It is easily seen that

[TABLE]

Since $\lambda_{K+1}(\Sigma)=0$ and $\lambda_{K}(\Sigma)\geq CM/p$ by Lemma 2, we obtain that

[TABLE]

It remains to bound $\|\hat{\Sigma}-\Sigma\|_{F}$. Note that we can decompose $\|\hat{\Sigma}-\Sigma\|_{F}$ as follows:

[TABLE]

where

[TABLE]

[TABLE]

and $\tilde{Z}$ is defined in Lemma 1. Lemmas 3 and 4 imply

[TABLE]

Since $Y_{m}$ and $Y_{m^{\prime}}$ are independent given $X$ and $Z$, $\mathbb{E}\epsilon^{\otimes 2}$ is a diagonal matrix. Combined with Assumption 3, we have

[TABLE]

The rest of the proof is to bound the last term $I_{4}$. By writing

[TABLE]

we obtain via Taylor expansion that

[TABLE]

where

[TABLE]

[TABLE]

is defined in the proof of Theorem 1 with $\bar{\phi}_{j}=\mathbb{E}(\phi_{j})$, and

[TABLE]

where $t_{j}$ is some intermediate value between $F_{j}^{*}X$ and $\Theta_{j}X+B_{j}Z$. In addition, we have

[TABLE]

where $\tilde{\zeta}_{j}$ is defined in the same way as $\zeta_{j}$ with $F_{j}^{*}X$ replaced by some intermediate value between $F_{j}^{*}X$ and $\hat{F}_{j}X$. Summarizing all the terms above, for a fixed observation $(i)$, looking at the column vector constructed by combining all $1\leq j\leq M$, we have $\hat{\epsilon}=\sum_{s=1}^{6}J_{s}$, where $J_{1}=\epsilon$ and $J_{s}$ is the column vector consisting of $J_{sj}$ for $s\geq 2$. Recall that we have I_{4}=\big{|}\big{|}\mathbb{E}_{n}\big{[}(J_{1}+...+J_{6})^{\otimes 2}\big{]}-\mathbb{E}_{n}\big{[}J_{1}^{\otimes 2}\big{]}-\mathbb{E}_{n}\big{[}J_{2}^{\otimes 2}\big{]}\big{|}\big{|}_{F}. Thus, we have

[TABLE]

For each term on the right hand side, we consider the diagonal and off-diagonal terms separately. Using the superscript $(i)$ to explicitly denote the $i$-th observation, note that each diagonal term in $\mathbb{E}_{n}J_{s}J_{t}^{T}$ will have form \big{|}\sum_{i=1}^{n}J_{sj}^{(i)}J_{tj}^{(i)}\big{|} and each off-diagonal term will have form \big{|}\sum_{i=1}^{n}J_{sj}^{(i)}J_{tk}^{(i)}\big{|} where $1\leq j\neq k\leq M$. Similarly, each diagonal term in $\mathbb{E}_{n}J_{s}^{\otimes 2}$ will have form \big{|}\sum_{i=1}^{n}J_{sj}^{(i)2}\big{|} while each off diagonal term will have form \big{|}\sum_{i=1}^{n}J_{sj}^{(i)}J_{sk}^{(i)}\big{|} where $1\leq j\neq k\leq M$. We start from the off-diagonal terms. For $j\neq k$ and $s=3$, we have $\mathbb{E}_{n}J_{sj}J_{sk}=\mathbb{E}_{n}\zeta_{j}\zeta_{k}B_{j}(Z-\tilde{Z})^{\otimes 2}B_{k}^{T}$. Since $\zeta_{j}\zeta_{k}$ has mean 0 conditioned on $X,Z$ and $\zeta_{j}\zeta_{k}$ is $1/2$-sub-exponential, Lemma 10 implies that conditioned on $X,Z$,

[TABLE]

provided $(\log M)^{3}/n=O(1)$. By Lemma 1, we know that $\|B_{m}(Z-\tilde{Z})\|_{\psi_{2}}\leq c/\sqrt{p}$, which implies $B_{j}(Z^{(i)}-\tilde{Z}^{(i)})^{\otimes 2}B_{k}^{T}$ is sub-exponential with norm of order $1/p$. By the tail bound for the maximum of sub-exponential r.v, we can show that

[TABLE]

For $s=4$,

[TABLE]

We use the same logic to bound $A_{j}:=(1+\zeta_{j})^{2}(\bar{\phi}_{j}\Sigma^{-1/2}X)^{2}$. Again, $(1+\zeta_{j})^{2}$ is $1/2$-sub-exponential and Lemma 10 implies that conditioned on $X,Z$,

[TABLE]

provided $(\log M)^{3}/n=O(1)$, where we use the fact that $\|\bar{\phi}_{j}\Sigma^{-1/2}X^{(i)}\|_{\psi_{2}}\lesssim 1/p$. In addition,

[TABLE]

By the sub-Gaussian property, $\mathbb{E}(\bar{\phi}_{j}\Sigma^{-1/2}X)^{2}\lesssim 1/p^{2}$. The Bernstein inequality implies

[TABLE]

Therefore,

[TABLE]

which implies

[TABLE]

Finally, we obtain

[TABLE]

For $s=5$, we have

[TABLE]

where $A_{5j}=\tilde{\eta}^{2}_{j}(F_{j}^{*}X-\Theta_{j}X-B_{j}Z)^{4}$. Following a similar argument,

[TABLE]

where we know from the proof of Theorem 1 that $\|F_{j}^{*}X-\Theta_{j}X-B_{j}Z\|_{\psi_{2}}\lesssim p^{-1/2}$. We can similarly show that

[TABLE]

and therefore under the assumption $(\log M)^{5}/n=O(1)$,

[TABLE]

Finally, we obtain

[TABLE]

For $s=6$,

[TABLE]

where $A_{6j}=(1+\tilde{\zeta}_{j})^{2}\{(\hat{F}_{j}-F_{j}^{*})X\}^{2}$. Due to sample splitting, given $X,Z$ and $\hat{F}$, the r.v. $(1+\tilde{\zeta}_{j})^{2}$ is $1/2$-sub-exponential and Lemma 10 implies that

[TABLE]

where the last line follows from Theorem 2. Together with

[TABLE]

we obtain that

[TABLE]

as we assume $p\sqrt{\frac{\log(p\vee M)}{n}}=o(1)$. As a result,

[TABLE]

Next, consider $s=1$ and $t=2$ in (50). Again, given $X$ and $Z$, Lemma 10 (or the Bernstein inequality in Lemma 9) implies

[TABLE]

Since $\|B_{m}(Z-\tilde{Z})\|_{\psi_{2}}\leq c/\sqrt{p}$, by the tail bound for the maximum of sub-Gaussian r.v, we can show that

[TABLE]

For $s=1$ and $t=3$,

[TABLE]

For $s=1$ and $t=4$, since $\epsilon_{j}(1+\zeta_{k})$ is mean 0 and $1/2$-sub-exponential given $X$ and $Z$, invoking Lemma 10 we have

[TABLE]

For $s=1$ and $t=5$,

[TABLE]

For $s=1$ and $t=6$, due to sample slitting, conditioned on $X,Z$ and $\hat{F}$, $\epsilon_{j}(1+\tilde{\zeta}_{k})$ is mean 0 and $1/2$-sub-exponential. Therefore

[TABLE]

where the last line follows from Theorem 2, the assumption that $X$ is bounded, and the fact that $||v||_{1}\leq\sqrt{p}||v||_{2}$ for any $p$-vector $v$. For $s=2$ and $t=3$,

[TABLE]

For $s=2$ and $t=4$,

[TABLE]

We notice that by the definition of $\tilde{Z}$, $B_{j}(Z-\tilde{Z})\bar{\phi}_{k}\Sigma^{-1/2}X$ is mean 0 and sub-exponential with sub-exponential norm of order $p^{-3/2}$. Given $X$ and $Z$, $\zeta_{k}$ is sub-exponential with bounded sub-exponential norm. Thus,

[TABLE]

and

[TABLE]

Combining above two terms,

[TABLE]

For $s=2$ and $t=5$,

[TABLE]

Since $\{B_{j}(Z-\tilde{Z})\}^{2}$ is sub-exponential with norm of order $1/p$, we have

[TABLE]

and $\mathbb{E}\{B_{j}(Z-\tilde{Z})\}^{2}\leq C/p$, which implies $\mathbb{E}_{n}\{B_{j}(Z-\tilde{Z})\}^{2}\lesssim 1/p$. The previous argument (for $s=5$) implies $\mathbb{E}_{n}\tilde{\eta}^{2}_{k}(F_{k}^{*}X-\Theta_{k}X-B_{k}Z)^{4}\lesssim 1/p^{2}$. As a result,

[TABLE]

Similarly, for $s=2$ and $t=6$,

[TABLE]

For $s=3$ and $t=4$,

[TABLE]

where we use the same argument used for $s=2$ and $t=4$ as $\zeta_{j}(1+\zeta_{k})$ is mean 0 and $1/2$-sub-exponential.

For $s=3$ and $t=5$, $\zeta_{j}\tilde{\eta}_{k}$ is mean 0 and $1/2$-sub-exponential given $X,Z$, and thus

[TABLE]

For $s=3$ and $t=6$, due to sample splitting, $\zeta_{j}(1+\tilde{\zeta}_{k})$ is mean 0 and $1/2$-sub-exponential given $X,Z$ and $\hat{F}$, and thus

[TABLE]

For $s=4$ and $t=5$,

[TABLE]

where the last step follows from the analysis for $s=4$ and $s=5$ above.

For $s=4$ and $t=6$, we apply the same argument to derive

[TABLE]

Finally, for $s=5$ and $t=6$

[TABLE]

We have considered all the off-diagonal terms on the right hand side of (50). For the diagonal terms, $j=k$, we have the following bounds. For brevity, we skip the intermediate steps.

[TABLE]

[TABLE]

It can be seen that the leading error among all the diagonal terms is

[TABLE]

By collecting the order of the errors in (51)–(69), the leading error among all the off-diagonal terms is

[TABLE]

Plugging these above results into (50), we derive

[TABLE]

since there are $M$ diagonal terms and $O(M^{2})$ off diagonal terms and thus we multiply the corresponding square roots when computing the Frobenius norm for $I_{4}$. From (48) and (49), after removing redundant terms, we finally obtain under $p\sqrt{p\log(p\vee M)/n}\cdot\log(n\vee M)=o(1)$ (assumed in Theorem 2) and $\log(n\vee M)=O(p^{3})$ the following:

[TABLE]

This completes the proof. ∎

B.4 Proof for Theorem 4

Proof.

For $\hat{\Theta}:=\hat{P}_{B}^{\perp}\hat{F}$, we have the following derivation:

[TABLE]

By Theorems 3, 2 and 1, we have

[TABLE]

∎

B.5 Proof for Theorem 5

Proof.

For any $u\in\mathbb{R}^{M}$ and $v\in\mathbb{R}^{p}$ with $\|u\|_{2}=1$ and $\|v\|_{2}=1$, we can write

[TABLE]

For the first term $I_{1}$, by (47) we can rewrite $I_{1}$ as

[TABLE]

where $\xi_{1}^{(i)},\xi_{2}^{(i)}\in\mathbb{R}^{M}$ with the $m$th entry being $\bar{\epsilon}^{(i)}_{m}v^{T}\Delta_{m1}X^{(i)}$ and $\bar{\epsilon}^{(i)}_{m}v^{T}\Delta_{m2}X^{(i)}$. Here,

[TABLE]

where the same notation in (47) is used here. In the following, we will first apply the Berry–Esseen theorem to the first term in (71). Let $O_{i}=u^{T}P_{B}^{\perp}h^{(i)}$. The Berry–Esseen theorem yields that

[TABLE]

where $s_{n}^{2}=\sum_{i=1}^{n}\mathbb{E}O_{i}^{2}$ and $\rho_{i}=\mathbb{E}|O_{i}|^{3}$. Since $\mathbb{E}O_{i}^{2}\geq C$, we have $s_{n}\geq\sqrt{Cn}$. By Hölder’s inequality, $\rho_{i}\leq\{\mathbb{E}O_{i}^{4}\}^{3/4}$. For $\mathbb{E}O_{i}^{4}$, we first note that

[TABLE]

where the last step follows from the factor model assumption, $\|u\|_{1}\leq R_{n}$ and the fact that each entry of $B$ is bounded by $C_{4}$ and therefore $\|B\|_{\infty,2}\lesssim\sqrt{K}$. Recall that $h^{(i)}_{m}=\bar{\epsilon}^{(i)}_{m}v^{T}G^{-1}_{m}X^{(i)}$, where $\bar{\epsilon}^{(i)}_{m}$ is sub-exponential, and that $v^{T}G^{-1}_{m}X^{(i)}=v^{T}\Sigma_{X}^{-1/2}\tilde{G}^{-1}_{m}\tilde{X}^{(i)}$ is sub-Gaussian with bounded sub-Gaussian norm. Also, $\tilde{X}^{(i)}=\Sigma_{X}^{-1/2}X^{(i)}$ is sub-Gaussian with bounded sub-Gaussian norm, $\tilde{G}_{m}=\mathbb{E}(1+\zeta_{m}^{(i)})\tilde{X}^{(i)}\tilde{X}^{(i)T}$ satisfies (45), and $\lambda_{\min}(\Sigma_{X})\geq 1$ (WLOG). Writing $u^{T}P_{B}^{\perp}h^{(i)}=u^{T}h^{(i)}-u^{T}P_{B}h^{(i)}$, we can show that

[TABLE]

Therefore, (72) implies

[TABLE]

For the second term in (71), we have

[TABLE]

By the proof of Theorem 2,

[TABLE]

and

[TABLE]

As a result, we can show that

[TABLE]

Following a similar argument,

[TABLE]

which implies that

[TABLE]

For the third term in (71), we have

[TABLE]

Combined with (74), we obtain the Berry–Esseen bound for $I_{1}$,

[TABLE]

It remains to bound $I_{2}$ and $I_{3}$. We note that by Lemma 5

[TABLE]

where $\eta$ is defined in (82) and derived in Remark 3. We then have

[TABLE]

As a result, the contribution to the Berry–Esseen bound from $I_{2}$ and $I_{3}$ is given by

[TABLE]

and

[TABLE]

To show (31), we decompose

[TABLE]

where we denote the above two terms as $J_{1}$ and $J_{2}$. Let

[TABLE]

By the Markov inequality we have

[TABLE]

In addition, we can further decompose $J_{2}$ as

[TABLE]

where $\tilde{h}^{(i)}=(\tilde{h}^{(i)}_{1},...,\tilde{h}^{(i)}_{M})^{T}$ with $\tilde{h}^{(i)}_{m}=\bar{\epsilon}^{(i)}_{m}v^{T}\hat{G}^{-1}_{m}X^{(i)}$. We denote the above two terms as $J_{21}$ and $J_{22}$, respectively. We first note that

[TABLE]

where we use (73). With a slight abuse of notation, let $\hat{H}_{m}=\frac{1}{n}\sum_{i=1}^{n}(1+\hat{\zeta}^{(i)}_{m})\tilde{X}^{(i)}\tilde{X}^{(i)T}$, $\tilde{H}_{m}=\frac{1}{n}\sum_{i=1}^{n}(1+\zeta^{(i)}_{m})\tilde{X}^{(i)}\tilde{X}^{(i)T}$ and $H_{m}=\mathbb{E}(1+\zeta^{(i)}_{m})\tilde{X}^{(i)}\tilde{X}^{(i)T}$, where $\tilde{X}^{(i)}=\Sigma_{X}^{-1/2}X^{(i)}$. Recall from the proof of Theorem 2 that

[TABLE]

and

[TABLE]

and recall that the smallest and largest eigenvalues of $H_{m}$ are lower and upper bounded by constants as shown in (44). As a result,

[TABLE]

To bound the above term, we first consider

[TABLE]

where the third line follows from $\max_{i}\|\bar{\epsilon}^{(i)}\|_{\infty}\lesssim\log(M\vee n)$, (77) and (78), and the last line is from

[TABLE]

Therefore,

[TABLE]

To bound $J_{22}$, we can decompose $J_{22}$ as

[TABLE]

We first note that $\max_{i}\|\bar{\epsilon}^{(i)}\|_{\infty}\lesssim\log(M\vee n)$, and $\max_{i}\|\hat{\epsilon}^{(i)}\|_{\infty}\lesssim\log(M\vee n)$ since we have

[TABLE]

Following a similar argument for $J_{21}$, we have

[TABLE]

where we use (76) and (79) in the final step. Finally, we further note that by Lemma 5, with $\eta$ defined therein, we have

[TABLE]

Combined with (76), $\|u^{T}\hat{P}_{B}^{\perp}\|_{1}\lesssim R_{n}\sqrt{K}+R_{n}\eta\sqrt{p}\lesssim R_{n}\sqrt{K}$, where indeed we have $\eta\sqrt{p}=o(1)$ under the assumption $\frac{p}{\sqrt{M}}=o(1)$ and $p\sqrt{\frac{\log(p\vee M)}{n}}=o(1)$. For the second term in the decomposition of $J_{22}$, we can obtain that

[TABLE]

By the decomposition of $J_{22}$, we finally derive the bound

[TABLE]

Combining the bound for $J_{1}$, $J_{21}$ and $J_{22}$, we prove that

[TABLE]

This completes the proof. ∎

B.6 Supplementary Lemmas

Lemma 2**.**

Under Assumptions 1 - 4, we have

[TABLE]

for some fixed constant $C>0$.

Proof.

[TABLE]

where $(\Sigma_{Z}^{-1}+A^{T}A)^{-1/2}$ is defined since $\Sigma_{Z}^{-1}+A^{T}A$ is symmetric positive definite. Since $\Sigma_{Z}^{-1}+A^{T}A$ has rank $K$, the above is equal to the smallest eigenvalue of its transpose as $XY$ shares the same non-zero eigenvalues as $YX$.

[TABLE]

where the last line follows from Assumption 4. ∎

Lemma 3**.**

Under Assumption 2 and 3, in the regime of $\log M<n$, we have the following bound:

[TABLE]

for some fixed constant $C>0$. Thus, we have

[TABLE]

Proof.

The $(m,m^{\prime})$-th element of \frac{1}{n}\bm{\epsilon}\bm{\epsilon}^{T}-\mathbb{E}\big{[}\frac{1}{n}\bm{\epsilon}\bm{\epsilon}^{T}\big{]} is \sum_{i=1}^{n}\big{\{}\epsilon_{m}^{(i)}\epsilon_{m^{\prime}}^{(i)}-E[\epsilon_{m}^{(i)}\epsilon_{m^{\prime}}^{(i)}]\big{\}}/n. When $m=m^{\prime}$, $\epsilon_{m}^{(i)2}-\mathbb{E}[\epsilon_{m}^{(i)2}]$ is a centered $\frac{1}{2}$-sub-exponential random variable (defined in Definition 1) with a norm bounded by $\sigma^{4}_{\epsilon,\max}$. By the concentration inequality in Lemma 10, we have

[TABLE]

The same logic holds when $m\neq m^{\prime}$, except now $\epsilon_{m}^{(i)}\epsilon_{m^{\prime}}^{(i)}$ is a centered random variable. The upper bound for the average is exactly the same. Using the union bound over all $M^{2}$ elements in the matrix, we get the conclusion. ∎

Lemma 4**.**

Under Assumptions 1 - 4 and the regime of $K<n$, we have the following bound:

[TABLE]

for some fixed constant $C>0$. Thus, we have

[TABLE]

Proof.

Note that for a $K\times 1$ centered sub-Gaussian random vector $V$, if we define $\hat{\Sigma}_{V}:=\sum_{i=1}^{n}V^{(i)}V^{(i)T}/n$ and $\Sigma_{V}:=\mathbb{E}[VV^{T}]$, then we have the following bound from Koltchinskii and Lounici (2017):

[TABLE]

If we consider $V=Z-\tilde{Z}$, we have

[TABLE]

For the last step, we follow the proof of Lemma 2 to obtain

[TABLE]

So, we have the conclusion:

[TABLE]

This completes the proof. ∎

Remark 3**.**

While Theorem 3 controls the Frobenius norm of $\hat{P}_{B}-P_{B}$, it is not sharp enough to bound $\max_{j}\|(\hat{P}_{B}-P_{B})e_{j}\|_{2}$ required in the proof of Theorem 5. To this end, we provide the necessary bound in Lemma 5 and provide the derivation for the relevant rate, $\eta$, in this remark. We first apply spectral decomposition

[TABLE]

where $\hat{V}\in\mathbb{R}^{M\times M}$ and $\hat{D}^{2}\in\mathbb{R}^{M\times M}$ are the corresponding eigenvectors and eigenvalues (in non-increasing order). We define the estimator

[TABLE]

where $\hat{D}_{K}\in\mathbb{R}^{K\times K}$ is the square root of the top $K$ eigenvalues and $\hat{V}_{K}\in\mathbb{R}^{M\times K}$ is the matrix of corresponding eigenvectors. Then $\frac{1}{M}\hat{\Sigma}\hat{V}_{K}=\hat{V}_{K}\hat{D}_{K}^{2}$ and $\hat{V}_{K}=\frac{1}{M}\hat{\Sigma}\hat{V}_{K}\hat{D}_{K}^{-2}$, which implies

[TABLE]

Note that

[TABLE]

Define $H=\frac{1}{M}\mathbb{E}_{n}(Z-\tilde{Z})^{\otimes 2}B^{T}\hat{V}_{K}\hat{D}_{K}^{-1}\sqrt{M}$. We have

[TABLE]

Recall that $\Sigma=B(\Sigma_{Z}^{-1}+A^{T}A)^{-1}B^{T}$. By the proof of Lemma 2, we know $C_{1}M/p\leq\lambda_{k}(\Sigma)\leq C_{2}M/p$. Since $\lambda_{k}(\hat{\Sigma}/M)=\lambda_{k}^{2}(\hat{D})$, by Weyl’s inequality and the proof of Theorem 3, under the assumption $K=o(n)$, $p=o(\sqrt{M})$, we have

[TABLE]

which implies $\lambda_{k}(\hat{D})\asymp p^{-1/2}$. We further have

[TABLE]

where the third step follows from the proof of Theorem 3 by just taking the sum of $M$ elements in one column of the $M\times M$ matrix in consideration rather than all $M^{2}$ elements. The last step is due to $\{\log(M\vee p)\}^{5}=O(n)$. Indeed, the above bound holds uniformly over $1\leq j\leq M$. For simplicity, we denote

[TABLE]

Lemma 5**.**

Under the same assumptions as in Theorem 3 and $K=o(n)$, $p=o(\sqrt{M})$, we have

[TABLE]

Proof.

Let $\tilde{B}=BH$. Since we can also write $P_{B}=\tilde{B}(\tilde{B}^{T}\tilde{B})^{-1}\tilde{B}^{T}$, it is apparent from the identity $A^{-1}-B^{-1}=-A^{-1}(A-B)B^{-1}$ that

[TABLE]

We denote the above three terms on the right hand side as $I_{1},I_{2}$, and $I_{3}$, respectively. For $I_{1}$,

[TABLE]

We note that $\|B^{T}e_{j}\|_{2}\leq C_{4}$ by Assumption 4, $\|(B^{T}B)^{-1}\|_{op}\lesssim M^{-1}$ and $\|\tilde{B}-\hat{B}\|_{F}\lesssim\eta\sqrt{M}$ by (81). In addition, by Weyl’s inequality and Lemma 4,

[TABLE]

which implies $\|H^{-1}\|_{\text{op}}\lesssim p^{1/2}$. Combining all these results, we obtain that

[TABLE]

For $I_{2}$, note that the smallest eigenvalue of $\hat{B}^{T}\hat{B}$ is lower bounded by $M/p$ since for $\tilde{B}$ we have $\lambda_{\min}(\tilde{B}^{T}\tilde{B})=\lambda_{\min}(H^{T}B^{T}BH)\geq\lambda_{\min}(B^{T}B)\lambda_{\min}(H^{T}H)\gtrsim\sqrt{M/p}\sqrt{M/p}=M/p$, and by Weyl’s inequality we have $\lambda_{\min}(\hat{B}^{T}\hat{B})\geq\lambda_{\min}(\tilde{B}^{T}\tilde{B})-||\hat{B}^{T}\hat{B}-\tilde{B}^{T}\tilde{B}||_{\text{op}}$, where the last term is upper bounded by \big{(}||\hat{B}||_{\text{op}}+||\tilde{B}||_{\text{op}}\big{)}||\hat{B}-\tilde{B}||_{F}=O(\sqrt{M/p})\cdot o(\sqrt{M/p})=o(M/p) in the assumed regime of $p=o(\sqrt{M})$. Combined with the following:

[TABLE]

we derive

[TABLE]

Finally, we can show that

[TABLE]

This completes the proof. ∎

Appendix C Additional Technical Results for Remark 1

In this section, we present some results on the eigenvalue ratio used in Remark 1. Recall that \Sigma=B\big{(}\Sigma_{Z}^{-1}+A^{T}A\big{)}^{-1}B^{T}. We know from Assumption 4 that \lambda_{j}\big{(}\Sigma\big{)}\asymp M/p for $1\leq j\leq K$, and \lambda_{j}\big{(}\Sigma\big{)}=0 for $K+1\leq j\leq M$. Similar to the derivation in the proof for Theorem 3, we know by Weyl’s inequality that for all $1\leq j\leq M$,

[TABLE]

Assume that

[TABLE]

Then $M/p$ dominates all of the terms in the above rate and we conclude that $\lambda_{j}(\hat{\Sigma})\asymp M/p$ for $1\leq j\leq K$, and $\lambda_{j}(\hat{\Sigma})=o(M/p)$ for $K+1\leq j\leq M$. This implies that $\lambda_{j}(\hat{\Sigma})/\lambda_{j+1}(\hat{\Sigma})\asymp 1$ for $1\leq j\leq K-1$, and $\lambda_{K}(\hat{\Sigma})/\lambda_{K+1}(\hat{\Sigma})\rightarrow\infty$. This implies $\hat{K}\geq K$. While the eigenvalue ratio approach in general cannot imply $\hat{K}=K$ with high probability, the numerical results in Bing et al. (2022) show that the performance of the PCA based estimator is robust even if $\hat{K}$ is above $K$, as long as it’s in a reasonable range.

Appendix D Other Useful Definitions and Inequalities

D.1 Davis-Kahan Theorem for Statisticians

Lemma 6**.**

Let $\Sigma,\hat{\Sigma}\in\mathbb{R}^{p\times p}$ be symmetric with eigenvalues $\lambda_{1}\geq...\geq\lambda_{p}$ and $\hat{\lambda}_{1}\geq...\geq\hat{\lambda}_{p}$. Fix $1\leq r\leq s\leq p$ and assume that $\min(\lambda_{r-1}-\lambda_{r},\leavevmode\nobreak\ \lambda_{s}-\lambda_{s+1})>0$ where we define $\lambda_{0}=\infty$ and $\lambda_{p+1}=-\infty$. Let $d=s-r+1$ and let $V=[v_{r},v_{r+1},...,v_{s}]\in\mathbb{R}^{p\times d}$ and $\hat{V}=[\hat{v}_{r},\hat{v}_{r+1},...,\hat{v}_{s}]\in\mathbb{R}^{p\times d}$ have orthonormal columns satisfying $\Sigma v_{j}=\lambda_{j}v_{j}$ and $\hat{\Sigma}\hat{v}_{j}=\hat{\lambda}_{j}\hat{v}_{j}$ for $j=r,r+1,...,s$. Then there exists an orthogonal matrix $\hat{O}\in\mathbb{R}^{d\times d}$ such that

[TABLE]

The full presentation of the theorem and its proof can be found in Yu et al. (2015).

D.2 Maximal Inequality for Sub-Gaussian Random Variables

Lemma 7**.**

Let $X_{1},...,X_{n}$ be $n$ random variables such that $X_{i}\sim$subGaussian$(||X||^{2}_{\psi_{2}})$. Then,

[TABLE]

D.3 Other Concentration Inequalities

Lemma 8**.**

$\mathrm{(Hoeffding\leavevmode\nobreak\ Inequality)}$* Let $X_{1},...,X_{n}$ be independent, mean-zero sub-Gaussian random variables, and let $a=(a_{1},...,a_{n})\in\mathbf{R}^{n}$. Then, for every $t\geq 0$, we have*

[TABLE]

where $V=\underset{1\leq i\leq n}{\max}||X_{i}||_{\psi_{2}}$.

Lemma 9**.**

$\mathrm{(Bernstein\leavevmode\nobreak\ Inequality)}$* Let $X_{1},...,X_{n}$ be independent, mean-zero sub-exponential random variables, and let $a=(a_{1},...,a_{n})\in\mathbf{R}^{n}$. Then, for every $t\geq 0$, we have*

[TABLE]

where $V=\underset{1\leq i\leq n}{\max}||X_{i}||_{\psi_{1}}$.

The proof for this lemma can be found in Vershynin (2018). Note that if you set T=\min\Big{(}\frac{t^{2}}{V^{2}\cdot||a||_{2}^{2}},\frac{t}{V\cdot||a||_{\max}}\Big{)}, this can be rewritten as:

[TABLE]

Definition 1**.**

$(\alpha-\mathrm{Sub}$-$\mathrm{exponential\leavevmode\nobreak\ Random\leavevmode\nobreak\ Variables)}$ We say that a random variable $X$ is $\alpha$-sub-exponential if there exists $K_{\alpha}>0$ such that $P(|X|\geq t)\leq 2\exp(-\frac{t^{\alpha}}{K_{\alpha}^{\alpha}})$ for all $t\geq 0$. We define the $\alpha$-sub-exponential norm as ||X||_{\psi_{\alpha}}:=\underset{p\geq 1}{\sup}\leavevmode\nobreak\ \frac{1}{p^{1/\alpha}}\cdot\big{\{}E(|X|^{p})\big{\}}^{1/p}.

Note that it follows that $E(|X|^{2})\leq 2^{\frac{2}{\alpha}}\cdot||X||_{\psi_{\alpha}}^{2}$. Thus, for a fixed $\alpha\in(0,1]$, a centered $\alpha$-sub-exponential random variable with a finite $\alpha$-sub-exponential norm has a bounded variance.

Lemma 10**.**

$(\mathrm{Concentration\leavevmode\nobreak\ for\leavevmode\nobreak\ }\alpha-\mathrm{Sub}$-$\mathrm{exponential\leavevmode\nobreak\ Random\leavevmode\nobreak\ Variables)}$ Let $X_{1},...,X_{n}$ be independent, mean-zero $\alpha$-sub-exponential random variables satisfying $||X_{i}||_{\psi_{\alpha}}\leq V$, for some $\alpha\in(0,1]$. Let $a\in\mathbf{R}^{n}$. For any $t>0$, we have

[TABLE]

The proof for Lemma 10 can be found in Götze et al. (2021). In a similar fashion to the lemma above, this can be rewritten as:

[TABLE]

In particular, by taking $T=\log M$ and plugging in $a=[1/n,...,1/n]$, we obtain

[TABLE]

provided $(\log M)^{(2-\alpha)/\alpha}/n=O(1)$.

Lemma 11**.**

$(\mathrm{Concentration\leavevmode\nobreak\ of\leavevmode\nobreak\ the\leavevmode\nobreak\ Euclidean\leavevmode\nobreak\ Norm\leavevmode\nobreak\ of\leavevmode\nobreak\ }\alpha-\mathrm{Sub}$-$\mathrm{exponential\leavevmode\nobreak\ Random\leavevmode\nobreak\ Variables)}$ Let $X_{1},...,X_{n}$ be independent, mean-zero $\alpha$-sub-exponential random variables satisfying $||X_{i}||_{\psi_{\alpha}}\leq V$, for some $\alpha\in(0,1]$ and $E(X^{2}_{i})\leq w_{i}^{2}$. For an $n\times n$ matrix $Q$, let $A=Q^{T}Q=(a_{ij})$. Then for any $t>0$, we have

[TABLE]

The proof for Lemma 11 can also be found in Götze et al. (2021).

Appendix E NHANES dataset

We now apply our G-hive procedure to the 2017-2018 NHANES dataset through the R package nhanesA (Schur, 2014). This is a public dataset that consists of a combination of demographic data (age, gender, etc.), examination data (height, weight, etc.), questionnaire responses (\sayDo you now smoke cigarettes every day, some days, or not at all?), and laboratory results (quantity of biomarkers taken from blood or urine samples). We take a subset of this data such that each response variable is binary and each covariate is continuous, and apply our data driven G-hive method on it. We focused on five binary responses corresponding to whether the individuals were diagnosed with depression, hypertension, arthritis, diabetes, and osteoporosis. We included age, income-to-poverty-ratio, systolic blood pressure (mmHg), body mass index (BMI, kg/$\text{m}^{2}$), and fasting glucose levels (mg/dL) as the observed covariates. After removing observations that have missing values for any of the response variables or covariates, we were left with $n=320$ valid observations. All of the observed covariates were standardized prior to the data analysis to prevent variables that are on a larger scale from dominating the model. The results with data driven G-hive for the coefficients $\widehat{\Theta}$ are presented in Table 3.

While it is infeasible in general to discern the ground truth coefficient values for real data analyses, it is apparent that many coefficients are well aligned with basic knowledge of the relationships between the covariates and the response variables. For instance, fasting glucose levels are known to be higher in individuals with diabetes (ElSayed et al. (2025)), and this is reflected in the large positive coefficient value of $0.9332$. Also, it is well known that BMI is a good predictor of diabetes (Hu et al. (2001)), and this is consistent with our relatively large positive value of $0.4284$. Another example of well-alignment with the medical literature is the large positive coefficient value of $0.4070$ that relates systolic blood pressure to hypertension as hypertension is normally diagnosed with a combination of systolic blood pressure values and diastolic blood pressure values (Parikh et al. (2008)). Lastly, age is shown to be positively correlated with arthritis, which is consistent with common knowledge in medicine as well (Elgaddal et al. (2024)).

E.1 Deconfounding benefits of G-HIVE on the NHANES dataset

We also present real data analysis results that highlight the ability of G-hive to account for model misspecification bias in the context of confounding. Recall that our method assumes that $Y$ depends on $X,Z$ through the GLM in (1) and that $X$ depends on the hidden variable $Z$ through the factor model $X=AZ+W$ in (2). This has a very natural connection to the basic confounded model depicted in Figure 3(a). Because $Z$ affects $X$ and $Y$, the observational association $P(Y|X)$ mixes the effect of $X\rightarrow Y$ and the spurious flow through $Z$ (Pearl (2009)). We utilize the same NHANES dataset from the previous section with slightly different variables that yield $n=230$ viable observations (with no missing values, etc.). We focus on two explanatory variables, \sayAge and \sayIncome to clearly see the effects of confounding from a hidden variable (\sayAge) and deconfounding with G-hive. Figure 3(b) shows the basic confounded model with \sayHypertension as an example response variable. Figure 3(b) is reasonable as it is well established that \sayAge has a positive effect on \sayIncome (Mincer (1974)) and that \sayAge has a positive effect on \sayHypertension (Parikh et al. (2008)). The same can be said about the effect of \sayAge on \sayDiabetes (Wilson et al. (2007)) and \sayAge on \sayArthritis (Elgaddal et al. (2024)), hence justifying similar figures with these response variables included instead. We run both naive mle and G-hive on the reduced model that just includes \sayIncome as the covariate, and we run naive mle on the full model that includes both \sayAge and \sayIncome as covariates. It is expected that the confounding will cause the coefficient corresponding to \sayIncome in the reduced model to appear more positive than it truly is in the full model. The coefficient values for each of these models are shown in Tables 4(b) and 5. We assume the latter model shows the \saytrue effect of \sayIncome on \sayHypertension, \sayDiabetes, and \sayArthritis after accounting for the effect of \sayAge. Comparing the coefficient values corresponding to \sayIncome on \sayArthritis between Table 4(a) and Table 4(b), it is apparent that G-hive’s $0.14$ is closer to the \saytrue value of $0.07$ compared to the more positively pushed naive mle value of $0.19$. Similarly, for \sayHypertension, G-hive’s $0.09$ is closer to the \saytrue value of $0.04$ compared to the positively shifted naive mle value of $0.16$. Lastly, for \sayDiabetes, G-hive’s $-0.38$ is closer to the \saytrue value of $-0.36$ compared to the positively shifted naive mle value of $-0.17$. This demonstrates that unlike naive mle, G-hive is able to account for confounding effects from the hidden variables and obtain estimates that are closer to the unconfounded effects even in real datasets.