High-Dimensional Overdispersed Generalized Factor Model with Application to Single-Cell Sequencing Data Analysis

In recent years, there has been a notable resurgence of high-dimensional factor models, which have proven to be valuable tools for analyzing complex datasets characterized by a large number of variables  [1, 2, 3]. These models have found widespread applications across various fields, including economics and finance for asset pricing [4], genomics for cell type identification  [5, 6], and social sciences for human ability assessment [7], among others. The versatility of high-dimensional factor models has positioned them as indispensable tools for addressing the challenges posed by intricate datasets and has paved the way for innovative research and analysis in diverse domains.

High-dimensional factor models provide a powerful framework for capturing the underlying structure and relationships within complex datasets. Through decomposing the observed variables into a reduced number of latent factors, these models enable dimension reduction and facilitate the extraction of meaningful information. The latent factors effectively capture the shared sources of variation across the variables, resulting in a more concise and interpretable representation of the data. In the current literature, high-dimensional factor models can be divided into two categories: linear factor models and nonlinear factor models.

Bai et al. [8] pioneered the high-dimensional linear factor model (LFM) and significantly advanced the field by establishing estimation theory and demonstrating the consistency of factor number selection. Since then, numerous studies have delved into high-dimensional LFMs [9, 1, 10, 3, 11]. LFMs exhibit excellent performance when the relationship between observed variables and factors is linear. However, for high-dimensional data with intricate dependencies, including nonlinearities, LFMs often fall short in terms of goodness of fit [2].

To address the limitations of high-dimensional LFMs, generalized factor models (GFMs) have been proposed as a class of models that utilize the exponential family of distributions to capture the nonlinear relationship between the high-dimensional observed variables and factors, such as Chen et al. [7], Wang [12] and Liu et al. [2]. Among these, Chen et al. [7] implicitly assumed a uniform exponential family distribution for all variables, which is not suitable for analyzing mixed-type data. In contrast, Liu et al. [2] and Wang [12] considered variable-specific distributions to model mixed-type data, where different variable types corresponded to different distributions. Unfortunately, existing nonlinear factor models are unable to account for overdispersion in mixed-type data, which may result in unsatisfactory estimation [13, 14]. Overdispersion is commonly encountered in practice, particularly in biomedical studies involving count responses, where the variability in the observed number of events often exceeds Poisson variability [13]. Additionally, overdispersion has been frequently observed in genomics, specifically in the analysis of single-cell RNA sequencing data [14, 15].

To overcome the limitations of existing models, we propose an overdispersed generalized factor model, called OverGFM, which is capable of simultaneously accounting for high-dimensional large-scale mixed data with overdispersion. Building upon the models proposed by Chen et al. [7] and Liu et al. [2], we formulate a hierarchical structure in OverGFM that incorporates an additional error term to explain overdispersion that cannot be captured by factors alone. However, OverGFM introduces significant computational challenges stemming from multiple factors. Firstly, it incorporates two high-dimensional latent random matrices, which contribute substantially to the computational complexity. Moreover, the model’s inherent nonlinearity adds an additional layer of complexity. To address these challenges, we introduce a variational EM (VEM) algorithm for implementing our model. The VEM algorithm combines Laplace and Taylor approximations, providing iterative explicit solutions for the complex variational parameters. Notably, our proposed VEM algorithm exhibits a high computational efficiency with linear complexity concerning sample size and variable dimension. We have theoretically proved the convergence of the proposed VEM algorithm. Furthermore, we develop a criterion based on the singular value ratio to determine the number of factors. In simulation studies, OverGFM showed improved estimation accuracy and remarkable computational efficiency in comparison to existing methods. Finally, we employed OverGFM to analyze two sets of single-cell sequencing data. The results unequivocally showcase its capacity in delivering invaluable biological insights within the genomics field, alongside its impressive computational scalability when addressing vast and intricate datasets.

The remaining sections of the paper are organized as follows. In Section 2, we provide an introduction to the model setup of OverGFM. Next, in Section 3, we present the estimation method, specifically focusing on the variational EM algorithm of OverGFM, as well as the procedure for selecting tuning parameter. To evaluate the performance of OverGFM, we conduct simulation studies in Section 4 and analyze real data in Section 5. In Section 6, we briefly discuss potential avenues for further research in this field. Technical proofs and additional numerical results are provided in the Supplementary Materials. Furthermore, we have seamlessly integrated OverGFM into an efficient and user-friendly R package, conveniently accessible at https://github.com/feiyoung/GFM.

Suppose that the observations $\{{\bf x}_{i},i=1,\cdots,n\}$, are independent and identically distributed (i.i.d.), where ${\bf x}_{i}=(x_{i1},\cdots,x_{ip})^{{}^{\mathrm{T}}}$ are variables of mixed types including continuous, binary, count variables, etc. Without loss of generality, we assume that there are $d$ variable types, and the index set of variables for each type $s$ is denoted by $G_{s},s=1,\cdots,d$. We consider an overdispersed generalized factor model given by a hierarchical formulation,

where $EF(\cdot)$ is an exponential family distribution and $g_{s}()$ is called mean function for variable type $s$. For example, if $x_{ij}$ is a continuous variable, then $EF(g_{s}(y_{ij}))=N(y_{ij},0)$, a degenerated normal distribution, i.e., $x_{ij}=y_{ij}$, and $g_{s}(y)=y$; if $x_{ij}$ is a count variable, then $EF(g_{s}(y_{ij}))=Poisson(\exp(y_{ij}))$ and $g_{s}(y)=\exp(y)$; and if $x_{ij}$ is a binary variable, then $EF(g_{s}(y_{ij}))=Bernoulli(\frac{1}{1+\exp(-y_{ij})})$ and $g_{s}(y)=\frac{1}{1+\exp(-y)}$. $a_{i}$ is a known offset term for unit $i$, ${\bf f}_{i}\in\mathbb{R}^{q}$ is a vector called latent factors, and ${\bf b}_{j}$ is the corresponding loading vector and $\mu_{j}$ is an intercept. The most significant difference between OverGFM and existing GFMs [7, 2] is that OverGFM can account for the extra variations in $x_{ij}$ not explained by factors. This is done with $\varepsilon_{ij}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,\lambda_{j})$, which considers these extra variations called overdispersion [13, 14]. Numerical findings demonstrate that this model design gives OverGFM a performance edge over existing GFMs.

Similar to Liu et al. [2], we mainly consider three variable types: continuous, count and binomial variables since they are popular in practice and the estimation procedure and the corresponding algorithm can be established similarly for other types belonging to the exponential family. Without loss of generality, let us assume types 1–3 corresponds to continuous, count and binomial variables, and denote $p_{s}=|G_{s}|$, where $|\cdot|$ is the cardinality of a set. The models of (2)–(2) for these variable types are explicitly written as

where $n_{j}$ is the number of trials for the $j$-th variable such that $j\in G_{3}$. If $n_{j}=1$ for all $j$, the binomial variable $x_{ij}$ reduces to the Bernoulli variable with success probability $p_{ij}$. Model (4) is unidentifiable due to the unobservability of ${\bf f}_{i}$ [2]. Let ${\bf B}=({\bf b}_{1},\cdots,{\bf b}_{p})^{{}^{\mathrm{T}}}$ be the loading matrix, ${\bf F}=({\bf f}_{1},\cdots,{\bf f}_{n})^{{}^{\mathrm{T}}}$ be the latent factor matrix, and ${\bf H}=({\bf h}_{1},\cdots,{\bf h}_{n})^{{}^{\mathrm{T}}}$ be the realization values of ${\bf F}$, i.e., factor score matrix. To make models computationally identifiable, we follow Bai et al. [9] and Liu et al. [2] to impose two conditions on the factor score matrix and loading matrix: (A1) $\frac{1}{n}\sum_{i=1}^{n}{\bf h}_{i}={\bf 0}$ and $\frac{1}{n}{\bf H}^{\mathrm{T}}{\bf H}={\bf I}_{q}$, where ${\bf I}_{q}$ is a $q$-by-$q$ identity matrix; and (A2) ${\bf B}^{\mathrm{T}}{\bf B}$ is diagonal with decreasing diagonal elements and the first nonzero element in each column of ${\bf B}$ is positive.

Let $p=\sum_{s=1}^{3}p_{s}$, ${\bf X}=(x_{ij},i=1,\cdots,n,j=1,\cdots,p)\in\mathbb{R}^{n\times p}$ and ${\bf Y}=(y_{ij},i=1,\cdots,n,j\in G_{2}\cup G_{3})\in\mathbb{R}^{n\times(p_{2}+p_{3})}$. Denote ${{\bm{\theta}}}=(\mu_{j},{\bf b}_{j},\lambda_{j},j=1,\cdots,{p},{\bf h}_{i},i=1,\cdots,n)$ that is the collection of unknown model parameters. The conditional log-likelihood (conditional on the latent factor matrix ${\bf F}$) of models (3)–(4) is derived as

by omitting the constant independent of parameters. There are significant computational challenges associated with the fact that ${\bf Y}$ is a large random matrix and ${\bf H}$ is a large unknown matrix. In existing GFMs [7, 2], only the factor score matrix ${\bf H}$ is unobservable. As a result, the computational challenges were addressed by treating the latent factors as ”parameters” to maximize the conditional log-likelihood [7, 2]. However, this approach is not applicable to our overdispersed GFM because of the additional unobservable large random matrix ${\bf Y}$. Therefore, we consider ${\bf Y}$ as latent variables handled by expectation-maximization (EM) algorithm, while ${\bf H}$ is regarded as a high-dimensional matrix parameter to be estimated directly. The EM algorithm [16] is a powerful and well-developed framework for handling models with latent variables, and involves the posterior distribution of the latent variables in a key step. However, in our model, computing the posterior distribution $P({\bf Y}|{\bf X},{\bf H})$ is extremely challenging due to the high dimensionality of both ${\bf Y}$ and $({\bf X},{\bf H})$, as well as the presence of nonlinear terms for Poisson and binomial variables in the log-likelihood (5).

To make the posterior distribution tractable, we utilize a mean field variational family, $q({\bf Y})$, to approximate $P({\bf Y}|{\bf X},{\bf H})$:

Let ${\bm{\gamma}}=(\tau_{ij},\sigma^{2}_{ij},i=1,\cdots,n,j\in G_{2}\cup G_{3})$ that is the collection of unknown variational parameters. In the proposed algorithm, $\bm{\gamma}$ is solved to seek an optimal approximation in the sense that KL divergence of $q({\bf Y})$ and $P({\bf Y}|{\bf X},{\bf H})$ is minimized.

Next, we derive the evidence lower bound (ELBO) function, which is given by

where $\mathrm{E}_{q(y_{ij})}F(y_{ij})$ is taking the expectation of $F(y_{ij})$ with respect to the random variable $y_{ij}\sim N(\tau_{ij},\sigma^{2}_{ij})$. In the following, we present a variational EM algorithm designed to implement the model.

In this section, we showcase the effectiveness of the proposed OverGFM through simulation studies involving 200 realizations. We compare OverGFM with various state-of-the-art methods from the current literature. They include

• (1)

  Generalized factor model [2] implemented in the R package GFM;

• (2)

  Multi-response reduced-rank regression model  [MRRR, 20] implemented in rrpack R package;

• (3)

  Principal component analysis for data with mix of qualitative and quantitative variables [21], implemented in the R package PCAmixdata;

• (4)

  Generalized PCA (GPCA) [22] implemented in the R package generalizedPCA;

• (5)

  High-dimensional LFM [8] implemented in the R package GFM;

• (6)

  Poisson PCA [23] implemented in the R package PoissonPCA;

• (7)

  PLNPCA [24] implemented in the R package PLNmodels;

Among the aforementioned methods, methods (1)-(3) are capable of handling mixed-type data. Method (4) can only analyze single-type data, including continuous, count, and categorical types. Method (5) is specifically designed for analyzing continuous variables, widely recognized as a benchmark with broad applications, notably in economics [25, 1] and genomics [26, 5], whereas methods (6) and (7) excel at analyzing count data.

We evaluate OverGFM in a total of eight scenarios. In scenarios 1-3, our main focus is comparing OverGFM with methods (1)-(3) and LFM, as LFM is widely utilized in practical applications [3, 6]. We generate data with a mixed-type of three variable types for these scenarios. In scenarios 4 and 5, we investigate the performance of the proposed SVR method in selecting the number of factors and the estimation performance under misselected $q$, respectively. In scenario 6, we investigate the computational efficiency of OverGFM by comparing it with other methods. In scenario 7, we focus on special cases where data is generated using a combination of two variable types or a single variable type, aiming to compare OverGFM with methods (4), (6) and (7). In scenario 8, we delve into the interconnection between OverGFM and GFM. To conserve space, the results (Table S1–S3, Figure S2) pertaining to scenarios 7–8 are deferred to Appendix B of the Supplementary Materials. In implementing the compared methods, we maintain the default settings and solely adjust the argument for the number of factors/principal components (PCs). To facilitate a fair comparison, we set the number of factors/PCs to the true value for all methods.

In scenarios 1–3, we generate data from models (2) and (2), i.e., $x_{ij}|y_{ij}\sim EF(g_{s}(y_{ij}))$, and $y_{ij}=a_{i}+{\bf b}_{j}^{{}^{\mathrm{T}}}{\bf h}_{i}+\mu_{j}+\varepsilon_{ij}$, and consider the mix of three different variable types: continuous, count and binary variables, i.e., $g_{1}(y)=y,g_{2}(y)=\exp(y)$ and $g_{3}(y)=1/(1+\exp(-y))$. Without loss of generality, we set the offset $a_{i}=0$ for all $i$’s. We set the number of variables of these three variable types to $\lfloor\frac{p}{3}\rfloor,\lfloor\frac{p}{3}\rfloor$ and $p-2\lfloor\frac{p}{3}\rfloor$, respectively. Next, we generate $\breve{{\bf B}}=(\breve{b}_{jk})\in\mathbb{R}^{p\times q}$ with $\breve{b}_{jk}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,1)$. Let $\breve{{\bf B}}_{s}$ be the submatrix of loading for variable type $s$. We generate $\bar{{\bf B}}_{s}=\rho_{s}\breve{{\bf B}}_{s}$, then construct $\bar{{\bf B}}=(\bar{{\bf B}}_{1}^{\mathrm{T}},\bar{{\bf B}}_{2}^{\mathrm{T}},\bar{{\bf B}}_{3}^{\mathrm{T}})^{\mathrm{T}}$, where $\rho_{s}$ controls the signal strength of each variable type. To obtain the singular value decomposition (SVD), we decompose $\bar{{\bf B}}=U_{2}\Lambda_{2}V_{2}^{{}^{\mathrm{T}}}$. Next, we define ${\bf B}_{0}=U_{2}\Lambda_{2}$. We then generate $\breve{{\bf h}}_{i}$ from $N({\bf 0}_{q},(0.5^{|i-j|})_{q\times q})$ and denote $\breve{{\bf H}}=(\breve{{\bf h}}_{1},\cdots,\breve{{\bf h}}_{n})^{{}^{\mathrm{T}}}$, perform column orthogonality for $\breve{{\bf H}}$ to obtain $\bar{\bf H}$, and set ${\bf H}_{0}=\bar{\bf H}V_{2}^{{}^{\mathrm{T}}}/\sqrt{n}$ such that ${\bf H}_{0}^{{}^{\mathrm{T}}}{\bf H}_{0}/n={\bf I}_{q}$. Note that ${\bf H}_{0}$ and ${\bf B}_{0}$ satisfy the identifiable conditions (A1)–(A2) given in Section 2. Finally, We generate $\mu_{j}=0.4z_{j},z_{j}\sim N(0,1)$, and $\varepsilon_{ij}\sim N(0,\sigma^{2})$. In all scenarios, we fix the number of latent factors ($q$) to 6.

We assess the estimation accuracy of the intercept-loading matrix $\Upsilon=(\bm{\mu},{\bf B})\in\mathbb{R}^{p\times(q+1)}$ and the factor matrix ${\bf H}$ by utilizing the commonly-used trace statistic [27] that measures the distance of the column space spanned by two matrices. For the factor score matrix, the trace statistic, denoted as $\mathrm{Tr}(\widehat{\bf H},{\bf H}_{0})$, is defined as $\mathrm{Tr}(\widehat{\bf H},{\bf H}_{0})=\frac{\mathrm{Tr}({\bf H}_{0}\widehat{\bf H}^{{}^{\mathrm{T}}}(\widehat{\bf H}^{{}^{\mathrm{T}}}\widehat{\bf H})^{-1}\widehat{\bf H}^{{}^{\mathrm{T}}}{\bf H}_{0})}{\mathrm{Tr}({\bf H}_{0}^{{}^{\mathrm{T}}}{\bf H}_{0})}$. The trace statistic yields a value between 0 and 1, with a higher value indicating greater accuracy in estimation.

In addition, we investigate the performance of the proposed method in comparison to its competitors when the overdispersion mechanism is incorrectly specified and there are highly heavy tails present. The results (Figure S1) show that OverGFM is not only flexible to the overdispersion but also robust to the heavy-tail data and model misspecification, making it a highly attractive and favorable choice in practical applications; see Appendix B.1 in Supplementary Materials.

In this section, we showcase the successful application of OverGFM in analyzing single-cell sequencing data within the genomics field. This includes the utilization of OverGFM on both a single-cell RNA sequencing (scRNA-seq) dataset and a single-cell multimodal sequencing dataset. The results demonstrate the effectiveness and versatility of OverGFM in handling diverse genomics data types.

We have introduced a novel statistical model named OverGFM, designed for the analysis of high-dimensional overdispersed mixed-type data. This model proves particularly valuable when dealing with scenarios where both the sample size (n) and variable dimension (p) are substantial. To address the computational challenges of large-scale data, we have developed a computationally efficient variational EM (VEM) algorithm. The VEM algorithm exhibits linear computational complexity with respect to sample size and variable dimension, ensuring its scalability. It offers straightforward implementation with explicit iterative solutions for all parameters and guarantees convergence, supported by theoretical proofs. Moreover, we developed a singular value ratio based method to determine the number of factors. In a series of numerical experiments, we have demonstrated that OverGFM surpasses existing methods, achieving superior estimation accuracy while reducing computation time. This makes OverGFM a compelling choice for analyzing extensive mixed-type datasets. Additionally, our application of OverGFM in the analysis of scRNA-seq and SNARE-seq data has proved its efficacy in unveiling the underlying structure of complex genomics data. We are confident that OverGFM holds the potential to enable essential discoveries not only in the field of genomics but also across other scientific domains. Its versatility and efficiency make it a promising tool for data analysis in various research areas.

A straightforward extension for OverGFM involves managing the high-dimensional mixed-type data with additional high-dimensional covariates. This extension would enable the model to explore the association between mixed-type variables and the extra covariates while also considering the presence of unobserved latent factors that cannot be explained by the covariates alone. We plan to pursue this direction in our future work, as it holds the promise of further enhancing the model’s capabilities and broadening its applicability to a wider range of real-world data analysis scenarios.