Sparse Bayesian Factor Models with Mass-Nonlocal Factor Scores

Let ${\bf Y}\in\mathbb{R}^{n\times p}$ be the $p$-dimensional data matrix where $n$ represents the number of observations and $p$ the number of observed variable. Assuming independence across observations and identical distribution, we can define a generic form of a latent factor model, for each unit $i\in\{1,\cdots,n\}$, as following:

where $\Lambda\in\mathbb{R}^{n\times p}$ is the factor loading matrix, $\eta_{i}\in\mathbb{R}^{k}$ is the latent factor score for the $i$-th observation, and $\epsilon_{i}\sim\mathcal{N}_{p}(0,\Sigma)$ represents the idiosyncratic error matrix, with $\Sigma=\mbox{diag}(\sigma_{1}^{2},\cdots,\sigma_{p}^{2})$.

Traditional factor models assume that the factor score are standard normally distributed, but this assumption often fails to capture sparsity and heterogeneity in the factor scores. To address this, we define a mixture distribution for $\eta$, assigning a nonzero probability to exact zeros using a Dirac distribution with mass in zero.

To ensure a flexible distribution that avoids overlap with the spike while maintaining tails similar to a normal distribution, we adopt a pMOM prior (Figure 1) as the slab component of the mixture [26, 27]. Hence, we define a mass-nonlocal prior, introduced by Shi et al. [42], for $\eta_{ih}$, for each unit $i$ and factor $h$, as follows:

where $\delta_{0}(\cdot)$ is a Dirac measure with mass at zero, and the probability density of the pMOM distribution [26, 27] is defined as follows:

where $\psi>0$ is the scale parameter.

The prior elicitation for the hyperparameters in the distribution involved in the (1) and (2) are are specified as follows:

This formulation allows us to introduce a latent variable $Z_{ih}$ for each unit $i\in\{1,\dots,n\}$ and factor $h\in\{1,\dots,H\}$, indicating whether the factor score belongs to the spike or the slab component. Specifically $Z_{ih}$ is defines as a binary random variable with the following distribution

where $\theta_{h}$ is the probability for the corresponding factor score $\eta_{ij}$ to have a pMOM probability distribution, such that

We adopt the multiplicative gamma process shrinkage (MGPS) prior introduced by Bhattacharya and Dunson [7] for the factor loading matrix $\Lambda$. However, possible alternative are the Bayesian cumulative shrinkage introduced by Legramanti et al. [30] or a generalization of Bhattacharya and Dunson [7]’s prior proposed by Schiavon et al. [39].

The MGPS prior [7] is defined for each entry of the factor loading matrix $\Lambda$ as follows:

where $\{\phi_{jh}\}_{h\geq 1,j\in\{1,\dots,p\}}$ are global shrinkage parameters follow gamma distribution with degree of freedom $\nu$, while $\{\tau_{h}\}_{h\geq 1}$ are global shrinkage parameters for the corresponding $h$-th column, obtained from the cumulative product of $\{\delta_{h}\}_{h\geq 1}$. Each $\delta_{h}$ follows a gamma distribution. Under specific settings, $\{\tau_{h}\}_{h\geq 1}$ is stochastically increasing process, favoring more shrinkage for higher-indexed factors, facilitating sparsity in $\Lambda$. We adopt the values for the hyperparameters suggested by [7, 19].

A straightforward Gibbs sampler can be proposed for posterior computation, levering the conjugate prior with the exception of pMOM distribution that requests a Metropolis-Hasting step.

Following the steps in the algorithm 2, in each iteration $r=1,\dots,R$, we use the observed data $y$ to update the following parameters and random variable: the factor loading $\Lambda$, denoting with $\lambda_{j}$ the $j$-th row, for $j\in\{1,\dots,p\}$, the factor score $\eta$, denoting with $\eta_{i}$ the $i$-th row, for $i\in\{1,\dots,n\}$, the diagonal matrix used for the MGPS prior (4) $D_{j}^{-1}=\mathbf{diag}(\phi_{j1}\tau_{1},\cdots,\phi_{jk}\tau_{k})$, the parameters $\{\theta_{j}\}_{j\in\{1,\dots,k\}}$, and the pMOM parameters $\{\psi_{j}\}_{j\in\{1,\dots,k\}}$. We indicated with $y^{(j)}=(y_{1j},\cdots,y_{nj})^{\intercal}$ for $j=1,\cdots,p$ the $n$ vector for the $j$-the outcome variable, and with $k$ the number of factors.

The steps for posterior sampling are as follows:

1. 1.

   The loading matrix entries $\{\lambda_{j}\}_{j\in\{1,\dots,p\}}$ are sampled from the following posterior distribution:

   $\displaystyle f(\lambda_{j}|\tau,\Lambda,D,\sigma_{y})\sim\mathcal{N}_{k}\Big{\{}\left(D_{j}^{-1}+\sigma_{j}^{-2}\eta^{\intercal}\eta\right)^{-1}\eta^{\intercal}\sigma_{j}^{-2}y^{(j)},\left(D_{j}^{-1}+\sigma_{j}^{-2}\eta^{\intercal}\eta\right)^{-1}\Big{\}}.$

2. 2.

   The the MGPS prior introduces two parameter. First, the local shrinkage parameter $\phi_{jh}$, with the following poster distribution, for $j\in\{1,\dots,p\}$ and $h\in\{1,\dots,k\}$:

   $\displaystyle f(\phi_{jh}|\nu,\Lambda,\tau)\sim\text{Ga}\left(\frac{\nu+1}{2},\frac{\nu+\tau_{h}\lambda_{jh}^{2}}{2}\right).$

3. 3.

   Second, the global shrinkage parameter $\delta_{h}$, with posterior distribution defined as follows:

   $\displaystyle f(\delta_{h}|a,\tau,\phi,\Lambda)\sim\text{Ga}\Big{\{}a_{h}+\frac{p}{2}\left(k-h+1\right),1+\frac{1}{2}\sum_{l=1}^{k}\tau_{l}^{(h)}\sum_{j=1}^{p}\phi_{jl}\lambda_{jl}^{2}\Big{\}},$

   where $\tau_{l}^{(h)}=\prod_{t=1,t\neq h}^{l}\delta_{t}$ for $h=1,\cdots,k$.

4. 4.

   The factor score $\eta_{ih}$, conditional to the latent variable $Z_{ih}$, are sampled from:

   $\displaystyle f(\eta_{ih}|Z_{ih})=\begin{cases}\pi(\eta_{ih}|c,d)&\mbox{ if }Z_{ih}=1,\\ 0&\mbox{ if }Z_{ih}=0;\end{cases}$

   where $\pi(\eta_{ih}|-)$ indicates the following distribution:

   	$$\displaystyle\pi(\eta_{ih}|c,d)\propto\exp\bigg{\{}-c\left(\eta_{ih}-\frac{d}{c}\right)^{2}\bigg{\}}\eta_{ih}^{2};$$
   	$$\displaystyle\mbox{with }c=\frac{1}{2\psi}+\sum_{j=1}^{p}\frac{1}{2\sigma_{j}^{2}}\lambda_{jh}^{2}\mbox{ and }d=\sum_{j=1}^{p}\frac{1}{2\sigma_{j}^{2}}\lambda_{jh}\left(y_{ij}-\sum_{l\neq h}^{k}\lambda_{jl}\eta_{il}\right).$$

   Due to the non-conjugacy of $\pi(\eta_{ih}|c,d)$, we embedded a Metroplis-Hastings algorithm, which is implemented as follows:
   

   Input: Probability density $\pi(\eta_{ih}|c,d)$, initial state $\eta_{ih}^{0}$

   Output: Posterior samples from $\pi(\eta_{ih}|c,d)$

   for $m=1$ to $M$ do

          Generate a random candidate $\eta_{ih}^{\ast}\sim\mathcal{N}(\mu=\eta_{ih}^{m-1},\sigma)$;

          Calculate acceptance probability $r=\exp\left(\log\pi(\eta_{ih}^{\ast}|c,d)-\log\pi(\eta_{ih}^{m-1}|c,d)\right)$;

          Accept or reject:

          $\alpha=\min(1,r),Z=\mbox{Bern}(\alpha)$;

          $\eta_{ih}^{m}=Z\eta_{ih}^{\ast}+(1-Z)\eta_{ih}^{m-1}$.

   end for

5. 5.

   Sample latent variable $Z_{ih}$, for $i\in\{1,\dots,n\}$ $h\in\{1,\dots,k\}$, from $\mbox{Bern}\left(f(Z_{ih}=0|-)\right)$, where:

   	$\displaystyle f(Z_{ih}=0|-)$	$\displaystyle=\frac{f(Z_{ih}=0)f(y_{i}|Z_{ih}=0,-)}{f(Z_{ih}=1)f(y_{i}|Z_{ih}=1,-)+f(Z_{ih}=0)f(y_{i}|Z_{ih}=0,-)}$
   		$\displaystyle=\frac{f(Z_{ih}=0)}{f(Z_{ih}=0)+f(Z_{ih}=1)T}$

   where $T=K\sqrt{2\pi}H^{-\frac{1}{2}}(H^{-1}+M^{2})$, $H=\frac{1}{\psi}+\lambda_{h}^{\intercal}\Sigma^{-1}\lambda_{h}$, $M=\frac{1}{H}(y_{i}-\Lambda_{(-h)}\eta_{i(-h)})\Sigma^{-1}\lambda_{h}$ and $K=2\pi\psi^{-\frac{3}{2}}\exp\{\frac{1}{2}HM^{2}\}$. The $\Lambda_{(-h)}$ represents $p\times(k-1)$ matrix with $h$th column dropped, and $\eta_{i(-h)}$ denotes $k-1$ vector with $\eta_{ih}$ entry deleted.

6. 6.

   Sample the probability parameter $\theta_{h}$, for each factor $h\in\{1,\dots,k\}$, from:

   $\displaystyle f(\theta_{h}|Z_{h},a_{1},b_{1})=\mbox{Beta}\left(\sum_{i}^{n}Z_{ih}+a_{1},n-\sum_{i}^{n}Z_{ih}+b_{1}\right);$

   where $a_{1}$ and $b_{1}$ are the hyperparameters.

7. 7.

   The residual variance $\sigma_{j}^{2}$, for $j\in\{1,\dots,p\}$, is drawn from the posterior distribution:

   $\displaystyle f(\sigma_{j}^{-2}|a_{\sigma},b_{\sigma},y,\Lambda,\eta)=\text{Ga}\Big{\{}a_{\sigma}+\frac{n}{2},b_{\sigma}+\frac{1}{2}\sum_{i=1}^{n}\left(y_{i}j-\lambda_{j}^{\intercal}\eta_{i}\right)^{2}\Big{\}}.$

A key challenge in factor analysis is the identifiability of the factor loading matrix $\Lambda$. For any orthogonal matrix $P$, $\Lambda^{\prime}=\Lambda P$ and $\eta_{i}^{\prime}=P^{\intercal}\eta_{i}$ satisfy equation 1 ($\Lambda\eta_{i}=\Lambda^{\prime}\eta_{i}^{\prime}$). This rotational invariance creates a non-identifiability issue, as multiple equivalent solutions exist for the same model. A common solution is to constrain the loading matrix to a lower triangular matrix, as described by Lopes and West [33]. These constraints can be restrictive and are not always practical. Moreover, each chain of the Gibbs sampler may converge to a different rotational solution, and averaging these solutions can lead to results that are not interpretable or meaningful. To address this challenge, Aßmann et al. [1] developed an orthogonal post-processing (OP) approach that resolves rotational ambiguity in unconstrained MCMC samples. This method generates samples without imposing rotational constraints, then aligns them by computing orthogonal transformations for each iteration. The procedure involves solving an optimization problem to align the sampled factor loading matrices and remove the effect of orthogonal rotations. When model parameters are not constrained, the Gibbs sampler produces orthogonally mixed MCMC samples. Each sample $\eta^{r}$ is associated with a potentially different orthogonal transformation $\big{\{}\mathbf{Q}^{(r)}\big{\}}_{r=1}^{R}$. The OP algorithm aligns these samples by minimizing the difference between the aligned samples and a reference loading matrix. Specifically, it solves the following optimization problem: Starting from a sequence of $R$ draws from the posterior distribution of $\big{\{}\eta^{(r)}\big{\}}_{r=1}^{R}$, the algorithm addressed this problem by estimating the loading matrices via the following constrained optimization:

where $L$ is the loss function:

The optimization is carried out by two steps while detailed derivation of this solution can be found at Schönemann [40]:

As Aßmann et al. [1] suggested, the last iteration of the Gibbs sampler can be chosen as the initial value for $\eta^{\ast}$. This iterative procedure aligns the loading matrices consistently across MCMC samples, ensuring interpretability and identifiability of the final results.

Convergence is assessed using the 2-norm to measure the distance between the current estimate and the final aligned loading matrix. Although the procedure is iterative, it typically requires only a single iteration to achieve alignment.

Determining the optimal number of factors is a critical challenge in factor analysis. The goal is to retain a small number of factors while capturing enough variance to preserve the structure of the data. Common approaches include threshold-based methods and information criteria.

Threshold-based methods, such as Kaiser’s criterion or scree plots, rely on subjective thresholds, such as retaining factors with eigenvalues greater than one in the Principal Component Analysis (PCA). Although easy to implement, these methods can be sensitive to data structure, leading to inconsistent results. Information criteria, such as the Bayesian Information Criterion (BIC), involve fitting models with varying numbers of factors and selecting the optimal model. However, these methods are computationally expensive, especially for high-dimensional data, as the process must be repeated for each potential number of factors. Furthermore, maximum likelihood estimates may not exist in such settings, reducing the reliability of these criteria.

Bhattacharya and Dunson [7] introduced an adaptive strategy to address the challenge of determining the optimal number of factors. This method begins with an overestimate of the number of factors and progressively shrinks redundant columns in the factor loading matrix. At each iteration, columns with elements near zero are pruned, while additional columns are introduced if needed. The effective number of factors stabilizes after a burn-in period, and the final estimate is determined using the median or mode of these stabilized values. This approach avoids the need for subjective thresholds and offers computational efficiency. However, simulations indicate that it may overestimate the number of factors in small sample sizes and low-dimensional settings.

In our approach, we have successfully leveraged the sparsity of the factor score matrix by specifying a mass-nonlocal prior. Building on this, instead of focusing on the loading matrix, we utilize the sparsity information inherent in the factor score matrix. Specifically, the parameter $\theta_{h}$—see definition of the model distribution (2) and the prior (3)—reflects the degree of sparsity for each factor $h$. Factors with higher sparsity are less likely to represent meaningful variation in the data and can be pruned from the candidate set of factors.

Unlike methods that apply a hard threshold on the proportion of non-zero entries in matrix to filter out auxiliary factors, our approach incorporates uncertainty of factor selection, which brings much flexibility. The whole process starts with a conservative guess, treating it as the true number of factors, then we implement our method for estimation. The simulation results demonstrate that, even if the specified number is incorrect, our model remains robust and produces highly accurate estimates. To determine the final number of factors, we examine the posterior mean of $\theta_{h}$, with respect to each factor $h$, which naturally reveal the probability that each factor is contributing or influential in capturing the variance structure of the data. Based on these probabilities, we can apply appropriate criteria or thresholds to ultimately determine the number of factors and truncate both the loading and the score matrices accordingly. This adaptive truncation balances robustness with flexibility, ensuring that only meaningful factors are retained while avoiding the pitfalls of overestimation.

To evaluate the ability to recover the factor loading matrix $\Lambda$ and the factor score matrix $\eta$ in the four simulated scenarios we use the RV coefficient [36]. The RV coefficient ranges from 0 to 1, where values closer to 1 indicate greater similarity between the estimated and true matrices.

The results, shown in Figure 2, demonstrate that our proposed model—referred to as Mass-nonlocal — effectively recovers the factors. In Scenario 1,2, and 4 the RV coefficients for $\Lambda\Lambda^{T}$ and $\eta\eta^{T}$ are consistently above $0.95$. In Scenario 3, where $p$ is twice the size of $n$, the RV coefficients are slightly lower, with median values close to $0.9$ for $\Lambda\Lambda^{T}$ and $0.85$ for $\eta\eta^{T}$, but the performance still exceeds that of the MGPS model.

The BFMAN model outperforms in all the scenarios, particularly in Scenarios 3 and 4, where its RV coefficients are often below 0.85. These results highlight the superior ability of the proposed model to recover the structure of the data under varying levels of sparsity and dimensionality.

The strong performance of our proposed model is attributed to its ability to effectively capture the sparsity of the factor score matrix. Figure 3 illustrates the distribution of the estimated probabilities of nonzero entries in the factor score matrix for each factor across the 50 replicates. The red lines in the figure indicate the true simulated values reported in Table 1. The estimated probabilities align well with the true values, as evidenced by their inclusion within the interquartile range of the boxplots in most cases.

Moreover, the BFMAN model demonstrates superior performance in selecting the number of factors. Table 3 presents the results using two different thresholds for BFMAN and MGPS. Specifically, the threshold of $0.17$ for $\theta$—i.e., the probability of observing non-zero entries in the factor scores for a specific factor—is chosen as it corresponds to the mean of the assumed prior Beta distribution with shape parameters equal to $1$ and $5$. Across all four scenarios, BFMAN consistently selects values closer to the true number of factors used in the data generation process, whereas MGPS tends to select notably higher values.

We estimated the proposed model with nutrients data, where the matrix $\bf Y$, defined in Section 2.1, is a matrix $12,972$. We define $k$—i.e., the number of factors used to run the model—equal to $6$. However, the posterior mean of the percentage of non-zero for factors is, in decreasing order, $(0.70,0.42,0.41,0.38,0.19,0.05)$, therefor we decided to consider just the first five factors and to remove the last one due to the $95\%$ of the corresponding factor scores are zeros.

Figure 4 reports the estimated factor load associated with nutrients. Inherently, according to the nutritional literature and with the results obtained for the same dataset by De Vito and Avalos-Pacheco [14], we can name the first factor as ’plant food’ due to the high-valued loadings in particolar on soluble and insoluble dietary fiber, vegetable protein, potassium, manganese, natural folate and copper. The second factor, with high values for animal protein, Zinc, Vitamin B12, B6, B2, and B3, Cholesterol, and short-chain saturated fatty acids (SCSFA) is nominated as ’animal product’. The ’seafood’ pattern is characterized by eicosapentaenoic acid (EPA), docosapentaenoic acid (DPA), docosahexaenoic acid (DHA), and arachidonic acid. Factor 4 can be defined as ’processed food’, due to the significant level of long-chain saturated fatty acids (LCSFA), long-chain monounsaturated fatty acids (LCSFA), linoleic and linolenic acids, and total trans fatty acids. Although factor 5, named ’dairy product’, has high-valued loadings on short- and medium-chain saturated fatty acids (SCSFA and MCSFA), calcium, retinol, and vitamin D.

The identified five dietary patterns are incorporated into a logistic regression model to investigate their association with cardiovascular diseases while controlling for gender. The results indicate that individuals with processed foods in their diet exhibit a significantly higher probability of developing hypercholesterolemia.