# Generalized infinite factorization models

Lorenzo Schiavon
Department of Statistical Sciences, University of Padova,
Via Cesare Battisti 241, 35121 Padova, Italy,
<EMAIL_ADDRESS>Antonio Canale
Department of Statistical Sciences, University of Padova,
Via Cesare Battisti 241, 35121 Padova, Italy,
<EMAIL_ADDRESS>and David B. Dunson
Department of Statistical Science, Duke University,
Durham, North Carolina 27708, U.S.A.
<EMAIL_ADDRESS>
###### Abstract

Factorization models express a statistical object of interest in terms of a
collection of simpler objects. For example, a matrix or tensor can be
expressed as a sum of rank-one components. However, in practice, it can be
challenging to infer the relative impact of the different components as well
as the number of components. A popular idea is to include infinitely many
components having impact decreasing with the component index. This article is
motivated by two limitations of existing methods: (1) lack of careful
consideration of the within component sparsity structure; and (2) no
accommodation for grouped variables and other non-exchangeable structures. We
propose a general class of infinite factorization models that address these
limitations. Theoretical support is provided, practical gains are shown in
simulation studies, and an ecology application focusing on modelling bird
species occurrence is discussed.

Adaptive Gibbs sampling; Bird species; Ecology; Factor analysis; High-
dimensional data; Increasing shrinkage; Structured shrinkage.

##  1 Introduction

Factorization models are used routinely to express matrices, tensors or other
statistical objects based on simple components. The likelihood for data $y$
under a general class of factorization models can be expressed as
$L(y;\Lambda,\Psi,\Sigma)$, with $\Lambda=\\{\Lambda_{h},h=1,\ldots,k\\}$ a
$p\times k$ matrix,
$\Lambda_{h}=(\lambda_{1h},\ldots,\lambda_{ph})^{\mathrm{\scriptscriptstyle
T}}$ the $h$th column vector of $\Lambda$, $\Psi$ and $\Sigma$ additional
parameters, and $k$ a positive integer. This class includes Gaussian linear
factor models (Roweis and Ghahramani, 1999), exponential family factor models
(Jun and Tao, 2013), Gaussian copula factor models (Murray et al., 2013),
latent factor linear mixed models (An et al., 2013), probabilistic matrix
factorization (Mnih and Salakhutdinov, 2008), underlying Gaussian factor
models for mixed scale data (Reich and Bandyopadhyay, 2010), and functional
data factor models (Montagna et al., 2012). A fundamental problem is how to
choose weights for the components and the number of components $k$. This
article proposes a general class of Bayesian methods to address this problem.

Although there is a rich literature, selection of $k$ is far from a solved
problem. In unsupervised settings, it is common to fit the model for different
choices of $k$ and then choose the value with the best goodness-of-fit
criteria. For likelihood models, the Bayesian information criteria is
particularly popular. It is also common to use an informal elbow rule,
selecting the smallest $k$ such that the criteria improves only a small amount
for $k+1$. In specific contexts, formal model selection methods have been
developed. For example, taking a Bayesian approach, one can choose a prior for
$k$ and attempt to approximate the posterior distribution of $k$ using Markov
chain Monte Carlo; see Lopes and West (2004) for linear factor models, Miller
and Harrison (2018) for mixture models and Yang et al. (2018) for matrix
factorization. Although such methods are conceptually appealing, computation
can be prohibitive outside of specialized settings.

Due to these challenges it has become popular to rely on over-fitted
factorization models, which include more than enough components with shrinkage
priors adaptively removing unnecessary ones by shrinking their coefficients
close to zero. Such approaches were proposed by Rousseau and Mengersen (2011)
for mixture models and Bhattacharya and Dunson (2011) for Gaussian linear
factor models. The latter approach specifically assumes an increasing
shrinkage prior on the columns of the factor loadings matrix $\Lambda$.
Legramanti et al. (2020) recently modified this approach using a spike and
slab structure (Mitchell and Beauchamp, 1988) that increases the mass on the
spike for later columns.

Although over-fitted factorizations are widely used, there are two key gaps in
the literature. The first is a careful development of the shrinkage properties
of increasing shrinkage priors (Durante, 2017). Outside of the factorization
context and mostly motivated by high-dimensional regression, there is a rich
literature recommending specific desirable properties for shrinkage priors.
These include high concentration at zero to favor shrinkage of small
coefficients and heavy tails to avoid over shrinking large coefficients.
Motivated by this thinking, popular shrinkage priors have been developed
including the Dirichlet-Laplace (Bhattacharya et al., 2015) and horseshoe
(Carvalho et al., 2010). Current increasing shrinkage priors, such as
Bhattacharya and Dunson (2011), were not designed to have the desirable
shrinkage properties of these priors. For this reason, ad hoc truncation and
use of horseshoe/Dirichlet-Laplace can outperform increasing shrinkage priors
in some contexts; for example, this was the case in Ferrari and Dunson (2020).

A second gap in the literature on over-fitted factorization priors is the lack
of structured shrinkage. The focus has been on priors for $\Lambda$ that are
exchangeable within columns, with the level of shrinkage increasing with the
column index. However, it is common in practice to have meta covariates
encoding features of the rows of $\Lambda$. For example, the rows may
correspond to different genes in genomic applications or species in ecology.
There is a rich literature on incorporating gene ontology in statistical
analyses of genomic data; refer, for example to Thomas et al. (2009). In
ecology it is common to include species traits in species distribution models
(Ovaskainen and Abrego, 2020). Beyond the Bayesian literature, it is common to
include structured penalties, with the grouped Lasso (Yuan and Lin, 2006) a
notable example.

Motivated by these deficiencies of current factorizations priors, this article
proposes a broad class of generalized infinite factorization priors, along
with corresponding theory and algorithms for routine Bayesian implementation.

##  2 Generalized infinite factor models

###  2.1 Model specification

Suppose that an $n\times p$ data matrix $y$ is available. In our motivating
application, $y_{ij}$ is a binary indicator of occurrence of bird species $j$
($j=1,\ldots,p$) in sample $i$ ($i=1,\ldots,n$). Consider the following
general class of models,

$\displaystyle y_{ij}=t_{j}(z_{ij}),\quad
z_{i}=\Lambda\eta_{i}+\epsilon_{i},\quad\epsilon_{i}\sim f_{\epsilon},$ (1)

with $\Lambda$ a $p\times k$ loadings matrix, $\eta_{i}$ a $k$ dimensional
factor with diagonal covariance matrix
$\Psi=\text{diag}(\psi_{11},\ldots,\psi_{kk})$, $\epsilon_{i}$ a
$p$-dimensional error term independent of $\eta_{i}$, and the function
$t_{j}:\Re\to\Re$, for $j=1\dots,p$. We refer to this class as generalized
factorization models. Class (1) includes most of the cases mentioned in
Section 1. When $\epsilon_{i}$ and $\eta_{i}$ are Gaussian random vectors and
$t_{j}$ is the identity function, model (1) is a Gaussian linear factor model.
With similar assumptions for $\epsilon_{i}$ and $\eta_{i}$, and assuming
$t_{j}=F_{j}^{-1}(\Phi(z_{ij}))$, with $\Phi(z_{ij})$ the Gaussian cumulative
distribution function, model (1) is a Gaussian copula factor model (Murray et
al., 2013). Exponential family factor models (Jun and Tao, 2013),
probabilistic matrix factorization (Mnih and Salakhutdinov, 2008) and
underlying Gaussian models for mixed scale data (Reich and Bandyopadhyay,
2010) can be obtained by appropriately defining the elements in (1), whereas
multivariate response regression models belong to this framework when
$\eta_{i}$ is known.

The matrix $\Omega=\text{var}(z_{i})$ can be expressed as
$\Omega=\Lambda\Psi\Lambda^{T}+\Sigma,$ where
$\Sigma=\text{var}(\epsilon_{i})$. Following common practice in Bayesian
factor analysis (Bhattacharya and Dunson, 2011), we avoid imposing
identifiability constraints on $\Lambda$ and assume $\Psi$ is pre-specified.
Our focus is on a new class of generalized infinite factor models induced
through a novel class of priors for $\Lambda$ that allows infinitely many
factors, $k=\infty$. In particular, we let

$\lambda_{jh}\mid\theta_{jh}\sim
N(0,\theta_{jh}),\quad\theta_{jh}=\tau_{0}\,\gamma_{h}\phi_{jh},\quad\tau_{0}\sim
f_{\tau_{0}},\quad\gamma_{h}\sim f_{\gamma_{h}},\quad\phi_{jh}\sim
f_{\phi_{j}},$ (2)

where $f_{\tau_{0}}$, $f_{\gamma_{h}}$, and $f_{\phi_{j}}$ are supported on
$[0,\infty)$ with positive probability mass on $(0,\infty)$. The local
$\phi_{jh}$, column-specific $\gamma_{h}$, and global $\tau_{0}$ scales are
all independent a priori. We let $N(0,0)$ denote a degenerate distribution
with all its mass at zero. Expression (2) induces a class of scale-mixture of
Gaussian shrinkage priors (Polson and Scott, 2010) for the loadings. Although
we allow infinitely many columns in $\Lambda$, (2) induces a prior for
$\Omega$ supported on the set of $p\times p$ positive semi-definite matrices
under mild conditions reported in Proposition S1 in the Supplementary
Material.

Differently from most of the existing literature on shrinkage priors, we want
to define a non-exchangeable structure that includes meta covariates $x$
informing the sparsity structure of $\Lambda$. In our context, meta covariates
provide information to distinguish the $p$ different variables as opposed to
traditional covariates that serve to distinguish the $n$ subjects. Letting $x$
denote a $p\times q$ matrix of such meta covariates, we choose $f_{\phi_{j}}$
not depending on the index $h$ and such that

$E(\phi_{jh}\mid\beta_{h})=g(x_{j}^{T}\beta_{h}),\quad\beta_{h}=(\beta_{1h},\ldots,\beta_{qh})^{\mathrm{\scriptscriptstyle
T}},\quad\beta_{mh}\sim f_{\beta}\quad(m=1,\ldots,q)$ (3)

where $g:\Re\to{\cal A}\subset\Re_{+}$ is a known smooth one-to-one
differentiable link function, $x_{j}=(x_{j1},\ldots,x_{jq})^{T}$ denotes the
$j$th row vector of $x$, and $\beta_{h}$ are coefficients controlling the
impact of the meta covariates on shrinkage of the factor loadings in the $h$th
column of $\Lambda$.

To illustrate the usefulness of (3), consider the previously introduced
ecological study and suppose
$x_{j}=\\{1,\mathbbm{1}{(\kappa_{j}=2)},\ldots,\mathbbm{1}{(\kappa_{j}=q)}\\}^{T}$,
where $\kappa_{j}\in\\{1,\ldots,q\\}$ denotes the phylogenetic order of
species $j$. Species of the same order may tend to have similarities that can
be expressed in terms of a shared pattern of high or low loadings on the same
latent factors.

Figure 1: Illustrative loadings matrix of an ecology application, where the
rows refer to ten bird species belonging to three phylogenetic orders. White
cells represent the elements of $\Lambda$ equal to zero, while blue and red
cells represent negative and positive values, respectively.

To illustrate this situation, we simulate a loadings matrix, displayed in Fig.
1, sampling from the prior introduced in Section 3 where
$\text{pr}(\lambda_{jh}=0)>\text{pr}(\phi_{jh}=0)>0$. The loadings within each
column are penalized basing on the group structure identified by the $q=3$
phylogenetic orders (Passeriformes, Charadriiformes, and Piciformes) of the
$p=10$ birds species considered. Our proposed prior allows for the possibility
of such structure while not imposing it. In the bird ecology application, $x$
can be defined to include not just phylogenetic placement of each bird species
but also species traits, such as size or diet (Tikhonov et al., 2020). Related
meta covariates are widely available, both in other ecology applications
(Miller et al., 2019) and in other fields such as genomics (Thomas et al.,
2009).

###  2.2 Properties

In this section we present some properties motivating the shrinkage process in
(2) and provide insight into prior elicitation. It is important to relate the
choice of hyperparameters to the signal-to-noise ratio, expressed as the
proportion of variance explained by the factors. Section S2.4 of the
Supplementary Material provides a study of the posterior distribution of the
proportion of variance explained; the posterior tends to be robust to
hyperparameter choice. Below we study key properties of our prior, including
an increasing shrinkage property, the ability of the induced marginal prior to
accommodate both sparse and large signals, and control of the multiplicity
problem in sparse settings. Proofs are included in the Appendix and in Section
S1 of Supplementary Material. This theory illuminates the role of
hyperparameters; specific recommendations of hyperparameter choice in practice
are illustrated under the model settings of Section 3.1.

To formalize the increasing shrinkage property, we introduce the following
definition.

######  Definition 1.

Letting $\Pi_{\Lambda}$ denote a shrinkage prior on $\Lambda$, $\Pi_{\Lambda}$
is a weakly increasing shrinkage prior if
$\text{var}(\lambda_{j(h-1)})>\text{var}(\lambda_{jh})$ for $j$ in
$1,\ldots,p$ and $h=2,\ldots,\infty$. $\Pi_{\Lambda}$ is a strongly increasing
shrinkage prior if $\text{var}(\lambda_{s(h-1)})>\text{var}(\lambda_{jh})$,
for $j,s$ in $\\{1,\ldots,p\\}$ and $h=2,\ldots,\infty$.

Weakly increasing shrinkage corresponds to the prior variance increasing
across columns within each row of $\Lambda$, while strongly increasing
shrinkage implies that the prior variance of any loading element is larger
than all elements with a higher column index. In the following Theorem, we
show that the process in (2) induces weakly increasing shrinkage under a
simple sufficient condition.

######  Theorem 1.

Expression (2) is a weakly increasing shrinkage prior under Definition 1 if
$\mbox{E}(\gamma_{h})>\mbox{E}(\gamma_{h+1})$ for any $h$.

Increasing shrinkage priors favor a decreasing contribution of higher indexed
columns of $\Lambda$ to the covariance $\Omega$. In addition to inducing a
flexible shrinkage structure that allows different factors to have a different
sparsity structure in their loadings, this allows one to accurately
approximate the likelihood $L(y;\Lambda,\Psi,\Sigma)$ by
$L(y;\Lambda_{H},\Psi_{H},\Sigma)$, with $\Lambda_{H}$ containing the first
$H$ columns of the infinite matrix $\Lambda$ and $\Psi_{H}$ the first $H$ rows
and columns of $\Psi$. To measure the induced truncation error of
$\Omega_{H}=\Lambda_{H}\Psi_{H}\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}}+\Sigma$, we use the trace of $\Omega$. The trace is justified by the fact
that the maximum error occurring in an element of $\Omega$ due to truncation
always lies along the diagonal and by the relation between difference of
traces and the nuclear norm, routinely used to approximate low rank
minimization problems (Liu and Vandenberghe, 2010). The following Proposition
provides conditions on prior (2) so that the under-estimation of $\Omega$ that
occurs by truncating decreases exponentially fast as $H$ increases.

######  Proposition 1.

Let $E(\tau_{0})$ and $E(\phi_{jh})$ be finite for $j=1,\ldots,p$ and
$h=1,\ldots,\infty$ and $E(\gamma_{h})=ab^{h-1}$ with $a>0$ and $b\in(0,1)$
for all $h=1,\ldots,\infty$. Let $c>0$ be a sufficiently large number such
that $c\geq\max_{h=1,\ldots,\infty}\psi_{hh}$. If

$m_{\Omega}=\min_{j=1,\ldots,p}\left[E(\sigma_{j}^{-2}),\,E\left\\{\left(\sum_{h=1}^{\infty}\psi_{hh}\lambda_{jh}^{2}\right)^{-1}\right\\}\right]<\infty,$

then for any $T\in(0,1)$,

$\text{pr}\left\\{\frac{\text{tr}(\Omega_{H})}{\text{tr}({\Omega})}\leq
T\right\\}\leq\bigg{(}\frac{1}{1-T}\bigg{)}\,ac\frac{b^{H}}{1-b}\,m_{\Omega}\,E(\tau_{0})\,\sum_{j=1}^{p}E(\phi_{j1}).$

The above increasing shrinkage properties can be satisfied by naive priors
that over-shrink the elements of $\Lambda$. It is important to avoid such
over-shrinkage and allow not only many elements that are $\approx 0$ but also
a small proportion of large coefficients. A similar motivation applies in the
literature on shrinkage priors in regression (Carvalho et al., 2010).
Borrowing from that literature, the marginal prior for $\lambda_{jh}$ should
be concentrated at zero to reduce mean square error by shrinking small
coefficients to zero but with heavy tails to avoid over-shrinking the signal.

To quantify the prior concentration of (2) in an $\epsilon$ neighbourhood of
zero, we can obtain

$\text{pr}(|\lambda_{jh}|>\epsilon)\leq\frac{E(\tau_{0})\,E(\gamma_{h})\,E(\phi_{jh})}{\epsilon^{2}}$
(4)

as a consequence of Markov’s inequality. Common practice in local-global
shrinkage priors chooses $\mbox{E}(\tau_{0})$ small while assigning a heavy-
tailed density to the local or column scales. In our case, (3) allows the
bound in (4) to be regulated by meta covariates $x$, while, under the
condition in Theorem 1, decreasing $E(\gamma_{h})$ with column index causes an
increasing concentration near zero, since $E(\phi_{jh})=E(\phi_{jl})$ for
every $h,l\in\\{1,\ldots,\infty\\}$. The means of the column and the local
scales control prior concentration near zero, while over-shrinkage can be
ameliorated by choosing $f_{\phi_{j}}$ or $f_{\gamma_{h}}$
($h=1,\ldots,\infty$) heavy tailed. The following Proposition provides a
condition on the prior to guarantee a heavy tailed marginal distribution for
$\lambda_{jh}$. A random variable has power law tails if its cumulative
distribution function $F$ has $1-F(t)\geq ct^{-\alpha}$ for constants $c>0$,
$\alpha>0$, and for any $t>L$ for $L$ sufficiently large.

######  Proposition 2.

If at least one scale parameter among $\tau_{0}$, $\gamma_{h}$ or $\phi_{jh}$
is characterized by a power law tail prior distribution, then the prior
marginal distribution of $\lambda_{jh}$ has power law tails.

An important consequence of the heavy tailed property is avoidance of over-
shrinkage of large signals. This is often formalized via a tail robustness
property (Carvalho et al., 2010). As an initial result, key to showing
sufficient conditions for a type of local tail robustness, we provide the
following Lemma on the derivative of the log prior in the limit as the value
of $\lambda_{jh}\to\infty.$

######  Lemma 1.

If at least one scale parameter among $\tau_{0}$, $\gamma_{h}$ or $\phi_{jh}$
has a prior with power law tails for any possible prior distribution of
$\beta_{h}$, then for any finite truncation level $H$,

$\lim_{\lambda\rightarrow\infty}\frac{\partial\log\\{f_{\lambda_{jh}\mid\Lambda_{-jh}}(\lambda)\\}}{\partial\lambda}=0$

where $f_{\lambda_{jh}\mid\Lambda_{-jh}}(\lambda)$ is the conditional
distribution of $\lambda_{jh}$ given the other elements of $\Lambda_{H}$.

The following definition introduces a type of local tail robustness property.

######  Definition 2.

Consider model (1) with factors $\eta$ known. Let $f_{\lambda_{jh}\mid
y,\eta,\Lambda_{-jh}}(\lambda)$ denote the posterior density of
$\lambda_{jh}$, given the data, conditional on any possible value of the other
elements of $\Lambda_{H}$ for any finite $H$, and let $\hat{\lambda}_{jh}$
denote the conditional maximum likelihood estimate of $\lambda_{jh}$ for any
possible value of the other elements of $\Lambda_{H}$. We say that the prior
on $\lambda_{jh}$ is tail robust if

$\lim_{\hat{\lambda}_{jh}\to\infty}\left|\hat{\lambda}_{jh}-\arg\max_{\lambda}f_{\lambda_{jh}\mid
y,\eta,\Lambda_{-jh}}(\lambda)\right|=0.$

For a given sample, $\hat{\lambda}_{jh}$ is a fixed quantity; the above limit
should be interpreted as what happens as the data support a larger and larger
maximum likelihood estimate. In order for tail robustness to hold, we need the
data to be sufficiently informative about the parameter $\lambda_{jh}$ and the
likelihood to be sufficiently regular; this is formalized as follows.

######  Assumption 1.

Let $L(y;\Lambda,\eta,\Sigma)$ denote the likelihood for data $y$
conditionally on latent variables $\eta$, let $l_{s}(\lambda)$ denote the
derivative function of the log-likelihood with respect to $\lambda_{jh}$, and
let $\mathcal{J}(\hat{\lambda}_{jh})$ denote the negative of the second
derivative of the log-likelihood with respect to $\lambda_{jh}$, evaluated at
the conditional maximum likelihood estimate $\hat{\lambda}_{jh}$. Then
$l_{s}(\lambda)$ is a continuous function for every $\lambda\in\Re$ and
$\mathcal{J}(\hat{\lambda}_{jh})\geq\nu(\hat{\lambda}_{jh})$, where
$\nu(\hat{\lambda}_{jh})$ is of order $O(1)$ as $\hat{\lambda}_{jh}\to\infty$.

This assumption can be verified for most of the cases mentioned in Section 1;
for example, for Gaussian linear factor models
$\mathcal{J}(\hat{\lambda}_{jh})$ is of order $O(1)$ with respect to
$\hat{\lambda}_{jh}$.

######  Theorem 2.

Under Assumption 1, if at least one scale parameter among $\tau_{0}$,
$\gamma_{h}$ or $\phi_{jh}$ is power law tail distributed for any possible
prior distribution of $\beta_{h}$, then the prior on $\lambda_{jh}$ is tail
robust under Definition 2.

As an additional desirable property, we would like to control for the
multiplicity problem within each column $\lambda_{h}$ of the loadings matrix,
corresponding to increasing numbers of false signals as dimension $p$
increases. This can be accomplished by imposing an asymptotically increasingly
sparse property on the prior, which is defined as follows.

######  Definition 3.

Let $|\text{supp}_{\epsilon}(\lambda_{h})|$ denote the cardinality of
$\text{supp}_{\epsilon}(\lambda_{h})=(j:|\lambda_{jh}|>\epsilon)$. Let
$s_{p}=o(p)$ such that $s_{p}\geq c_{s}\log(p)/p$ for some constant $c_{s}>0$.
We say that the prior on $\Lambda$ defined in (2) is an asymptotically
increasingly sparse prior if

$\lim_{p\to\infty}\text{pr}\\{|\text{supp}_{\epsilon}(\lambda_{h})|>a\,s_{p}\mid\gamma_{h},\tau_{0}\\}=0,\qquad\mbox{for
some constant $a>0$}.$

The quantity $|\text{supp}_{\epsilon}(\lambda_{h})|$ represents an approximate
measure of model size for continuous shrinkage priors and, conditionally on
$\beta_{h}$, $\gamma_{h}$, and $\tau_{0}$, it is a priori distributed as a sum
of independent Bernoulli random variables $\text{Ber}(\zeta_{\epsilon jh})$,
where

$\zeta_{\epsilon
jh}=\text{pr}(|\lambda_{jh}|>\epsilon\mid\beta_{h},\gamma_{h},\tau_{0})\leq\frac{\tau_{0}\,\gamma_{h}\,g(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})}{\epsilon^{2}}.$

We now provide sufficient conditions for an asymptotically increasingly sparse
prior, allowing regulation of the sparsity behaviour of the prior of the
columns of $\Lambda$ for increasing dimension $p$.

######  Theorem 3.

Consider prior (2) with $\phi_{jh}$ ($j=1,\ldots,p$) a priori independent
given $\beta_{h}$. If $\text{pr}\\{g(x_{j}^{{\mathrm{\scriptscriptstyle
T}}}\beta_{h})\leq\nu_{j}(p)\\}=1$, with $\nu_{j}(p)=O\\{\log(p)/p\\}$,
($j=1,\ldots,p)$, then the prior on $\Lambda$ is asymptotically increasingly
sparse under Definition 3.

The condition of the theorem is easily satisfied, for example, if $g$ is the
multiplication of a bounded function and a suitable offset depending on $p$ as
assumed in Section 3.1. The multiplicative gamma process (Bhattacharya and
Dunson, 2011) and cumulative shrinkage process (Legramanti et al., 2020) do
not satisfy the sufficient conditions of Theorem 3, and, furthermore, the
following lemma holds.

######  Lemma 2.

The multiplicative gamma process prior (Bhattacharya and Dunson, 2011) and the
cumulative shrinkage process prior (Legramanti et al., 2020) are not
asymptotically increasing sparse under Definition 3.

Although this Section has focused on properties of the prior, we find
empirically that these properties tend to carry over to the posterior, as will
be illustrated in the subsequent sections. For example, the posterior exhibits
asymptotic increasing sparsity; refer to Table 2 of Section 4, which shows
results for a novel process in our proposed class that is much more effective
than current approaches in identifying the true sparsity structure,
particularly when $p$ is large.

##  3 Structured increasing shrinkage process

###  3.1 Model specification

In this section we propose a structured increasing shrinkage process prior for
generalized infinite factor models satisfying all the sufficient conditions in
Propositions 1–2 and Theorems 2–3. Let $\text{Ga}(a,b)$ denote the gamma
distribution with mean $a/b$ and variance $a/b^{2}$. Following the notation of
Section 2.1, we specify

$\displaystyle\tau_{0}=1,\quad\gamma_{h}=\vartheta_{h}\rho_{h},\quad\phi_{jh}\mid\beta_{h}\sim\text{Ber}\\{\text{logit}^{-1}(x_{j}^{{\mathrm{\scriptscriptstyle
T}}}\beta_{h})\,c_{p}\\},$ (5)
$\displaystyle\vartheta_{h}^{-1}\sim\text{Ga}(a_{\theta},b_{\theta}),\quad
a_{\theta}>1,\quad\rho_{h}=\text{Ber}\left(1-\pi_{h}\right),\quad\beta_{h}\sim
N_{q}(0,\sigma_{\beta}^{2}I_{q}),$

where we assume the link $g(x)=\text{logit}^{-1}(x)c_{p}$, with
$\text{logit}^{-1}(x)=e^{x}/(1+e^{x})$ and $c_{p}\in(0,1)$ a possible offset.
The parameter $\pi_{h}=\mbox{pr}(\gamma_{h}=0)$ follows a stick-breaking
construction,

$\pi_{h}=\sum_{l=1}^{h}w_{l},\quad w_{l}=v_{l}\prod_{m=1}^{l-1}(1-v_{m}),\quad
v_{m}\sim\text{Be}(1,\alpha),$

with $\text{Be}(a,b)$ the beta distribution with mean $a/(a+b)$, such that
$\pi_{h+1}>\pi_{h}$ is guaranteed for any $h=1,\ldots,\infty$ and
$\lim_{h\rightarrow\infty}\pi_{h}=1$ almost surely. The prior expected number
of non degenerate $\Lambda$ columns is $E(\sum_{h=1}^{\infty}\rho_{h})=\alpha$
(Legramanti et al., 2020), suggesting setting $\alpha$ equal to the expected
number of active factors. The prior specification is completed assuming
$\Sigma=\text{diag}(\sigma_{1}^{2},\ldots,\sigma_{p}^{2})$ with
$\sigma_{j}^{-2}\sim\text{Ga}(a_{\sigma},b_{\sigma})$ for $j=1,\ldots,p$,
consistently with the literature. The hyperparameters can be chosen based on
one’s prior expectation of the signal-to-noise ratio, as $\sigma_{j}^{2}$ is
the contribution of the noise component to the total variance of the $j$th
variable. A sensitivity study in Section S2.4 of the Supplementary Material,
however, shows that posterior distributions tend to be robust to the
specification of $a_{\sigma},b_{\sigma}$. Regarding prior elicitation, we
recommend setting $b_{\theta}\geq a_{\theta}$ to induce a high enough
proportion of variance explained by the factor model. In Section S2.4 in the
Supplementary Materials we report empirical evidence of the effect of
different prior parameters on this quantity.

The above specification respects (2) and, consequently, the following
corollary holds.

######  Corollary 1.

The structured increasing shrinkage process defined in (5)

  * $i.$

is a strongly increasing shrinkage prior under Definition 1;

  * $ii.$

for any $T\in(0,1)$,

$\text{pr}\left\\{\frac{\text{tr}(\Omega_{H})}{\text{tr}({\Omega})}\leq
T\right\\}\leq\bigg{(}\frac{1}{1-T}\bigg{)}\,\frac{b^{H}}{1-b}\,\theta_{0}\,\frac{a_{\sigma}}{b_{\sigma}}\,\sum_{j=1}^{p}E(\phi_{j1}),$

with $b=\\{\alpha(1+\alpha)\\}^{-1}$ and $\theta_{0}=E(\vartheta_{h})$.

We conducted a simulation study on the posterior distribution of
$\\{\text{tr}(\Omega_{H})/\text{tr}({\Omega})\leq T\\}$ for varying
hyperparameters, and found that the results, reported in Section S2.4 of the
Supplementary Material, were quite consistent with our prior truncation error
bounds.

The prior concentration of the structured increasing shrinkage process in (5)
follows from (4):

$\text{pr}(|\lambda_{jh}|>\epsilon)\leq\frac{E(\vartheta_{h})\\{1-E(\pi_{h})\\}E(\phi_{jh})}{\epsilon^{2}}=\frac{\theta_{0}\,\\{\alpha/(1+\alpha)\\}^{h}}{\epsilon^{2}}\frac{c_{p}}{2}.$

In addition, the inverse gamma prior on $\vartheta_{h}$ implies a power law
tail distribution on $\gamma_{h}$ inducing robustness properties on
$\lambda_{jh}$ as formalized by the next corollary of Proposition 2 and
Theorem 2.

######  Corollary 2.

Under the structured increasing shrinkage process defined in (5)

  * $i.$

the marginal prior distribution on $\lambda_{jh}$ ($j=1,\ldots,p$;
$h=1,2,\ldots$) has power law tails;

  * $ii.$

under Assumption 1, the prior on $\lambda_{jh}$ ($j=1,\ldots,p$;
$h=1,2,\ldots$) is tail robust under Definition 2.

Finally, it is important to assess the joint sparsity properties of the prior
on each column of $\Lambda$. This is formalized in the following corollary of
Theorem 3.

######  Corollary 3.

If $c_{p}=O\\{\log(p)/p\\}$ the structured increasing shrinkage process
defined in (5) is asymptotically increasingly sparse under Definition 3.

###  3.2 Posterior computations

Posterior inference is conducted via Markov chain Monte Carlo sampling.
Following common practice in infinite factor models (Bhattacharya and Dunson,
2011; Legramanti et al., 2020; Schiavon and Canale, 2020) we use an adaptive
Gibbs algorithm, which attempts to infer the best truncation level $H$ while
drawing from the posterior distribution of the parameters. The value of $H$ is
adapted only at some Gibbs iterations by discarding redundant factors and, if
no redundant factors are identified, by adding a new factor by sampling its
parameters from the prior distribution. Convergence of the Markov chain is
guaranteed by satisfying the diminishing adaptation condition in Theorem 5 of
Roberts and Rosenthal (2007), by specifying the probability of occurrence of
an adaptive iteration $t$ as equal to $p(t)=\exp(\alpha_{0}+\alpha_{1}t)$,
where $\alpha_{0}$ and $\alpha_{1}$ are negative constants, such that
frequency of adaptation decreases.

The decomposition of $\gamma_{h}$ into two parameters $\rho_{h}$ and
$\vartheta_{h}$ allows one to identify the inactive columns of $\Lambda$,
corresponding to the redundant factors, as those with $\rho_{h}=0$, while
$H_{a}$ indicates the number of active columns of $\Lambda$. Consequently, at
the adaptive iteration $t+1$, the truncation level $H$ is set to
$H^{(t+1)}=H^{(t)}_{a}+1$ if $H^{(t)}_{a}<H^{(t)}-1$, and
$H^{(t+1)}=H^{(t)}+1$ otherwise. Given $H^{(t+1)}$, the number of factors of
the truncated model at iteration $t+1$, the sampler draws the model parameters
from the corresponding posterior full conditional distributions. The detailed
steps of the adaptive Gibbs sampler for the structured increasing shrinkage
prior in case of Gaussian or binary data are reported in the Supplementary
Material, as well as trace plots of the posterior samples for some parameters
of the model in Section 5 (see Section S3.2), showing good mixing.

###  3.3 Identifiability and posterior summaries

Non-identifiability of the latent structure creates problems in interpretation
of the results from Markov chain Monte Carlo samples. Indeed, both $\Lambda$
and $\eta$ are only identifiable up to an arbitrary rotation $P$ with
$PP^{T}=I_{k}$. This is a well known problem in Bayesian factor models and
there is a rich literature proposing post-processing algorithms that align
posterior samples $\Lambda^{(t)}$, so that one can then obtain interpretable
posterior summaries. Refer to McParland et al. (2014), Aßmann et al. (2016),
and Roy et al. (2019) for alternative post-processing algorithms in related
contexts.

Unfortunately, such post-hoc alignment algorithms destroy the structure we
have carefully imposed on the loadings in terms of sparsity and dependence on
meta covariates. Therefore, we propose a different solution to obtain a point
estimate of $\Lambda$ based on finding a representative Monte Carlo draw
$\Lambda^{(t)}$ consistently with the proposals of Dahl (2006) and Wade et al.
(2018) in the context of Bayesian model-based clustering. Specifically, we
summarize $\Lambda$ and $\beta=(\beta_{1},\beta_{2},\dots)$ through
$\Lambda^{(t^{*})}$ and $\beta^{(t^{*})}$ sampled at iteration $t^{*}$,
characterized by the highest marginal posterior density function
$f(\Lambda,\beta,\Sigma\mid y)$ obtained by integrating out the scale
parameters $\tau_{0},\gamma_{h},\phi_{jh}$ ($j=1,\ldots,p$, $h=1,\ldots$) and
the latent factors $\eta_{i}$ ($i=1,\ldots,n$) from the posterior density
function. Formally, we select the iteration $t^{*}\in\\{1,\ldots,T\\}$ such
that

$f(\Lambda^{(t^{*})},\beta^{(t^{*})},\Sigma^{(t^{*})}\mid
y)>f(\Lambda^{(t)},\beta^{(t)},\Sigma^{(t)}\mid y)\quad(t=1,\ldots,T),$

where $t=1,\ldots,T$ indexes the posterior samples. Under the structured
increasing shrinkage prior described in Section 3.1, these computations are
straightforward. The matrices
$\Lambda^{(t^{*})},\beta^{(t^{*})},\Sigma^{(t^{*})}$ are Monte Carlo
approximations of the maximum a posteriori estimator, which corresponds to the
Bayes estimator under $L_{\infty}$ loss. Although one can argue that
$L_{\infty}$ is not an ideal choice of loss philosophically in continuous
parameter problems, it nonetheless is an appealing pragmatic choice in our
context and is broadly used in other sparse estimation contexts, as in the
algorithm proposed by Ročková and George (2016) that similarly aims to recover
a strongly sparse posterior mode of an over-parameterized factor model.

##  4 Simulation experiments

We assess the performance of our structured increasing shrinkage prior
compared with current approaches (Bhattacharya and Dunson, 2011; Ročková and
George, 2016; Legramanti et al., 2020) through a simulation study. We have a
particular interest in inferring sparse and interpretable loadings matrices
$\Lambda$, but also assess performance in estimating the induced covariance
matrix $\Omega$ and number of factors. We generate synthetic data from four
scenarios based on different loadings structures. For each scenario we
simulate $R=25$ data sets with $n=250$ observations from $y_{i}\sim
N_{p}(0,\Lambda_{0}\Lambda_{0}^{T}+I_{p})\ (i=1,\ldots,n)$. In Scenario a, we
assume non sparse $\Lambda_{0}$, sampling the loadings $\lambda_{jh}$ from a
Gaussian distribution with mean zero, variance equal to
$\sigma^{2}_{\lambda}=1$ and ordering them to obtain decreasing variance over
the columns. To ensure that each element $\lambda_{jh}$ represents a signal,
we shifted them away from zero by $\sigma^{2}_{\lambda}/3$. In Scenario b we
remove the decreasing behaviour and introduce a random sparsity pattern
characterized by an increasing number of zero entries over the column index.
The loadings matrix for Scenario c is characterized by both the decreasing
behaviour over the columns of Scenario a and the random sparsity structure of
Scenario b. Finally, in Scenario d, while the decreasing behaviour is kept, we
induce a sparsity pattern dependent on a categorical and two continuous meta
covariates $x_{0}$. Details are reported in Section S2.2 of the Supplementary
Material.

For each scenario we consider four combinations of dimension and sparsity
level of $\Lambda_{0}$. We let $(p,k,s)\in\\{$(16,4,0.6), (32,8,0.4), (64,12,
0.3), (128,16,0.2)$\\}$, where $s$ is the proportion of non-zero entries of
$\Lambda$, with the exception of Scenario a where $s=1$. In these settings the
algorithm takes from $0.07$ to $0.73$ seconds of computational time per
iteration depending on the dimension $p$ and considering an R implementation
on an Intel Core i5-6200U CPU laptop with 15.8 GB of RAM. To estimate the
structured increasing shrinkage model, we set $x$ equal to the $p$-variate
column vector of $1$s, $\sigma_{\beta}=1$ and, consistently with Corollary 3,
$c_{p}=2e\log(p)/p$. In Scenario d we also estimate and compare a correctly
specified structured increasing model with $x=x_{0}$. For the method proposed
by Ročková and George (2016), we set the hyperparameters as suggested by the
authors, while for the remaining approaches, we follow the hyperparameter
specification and factor selection guidelines in Section 4 of Schiavon and
Canale (2020).

Table 1: Median and interquartile range of LPML and $E(H_{a}\mid y)$ in 25
replications of Scenario a for different combinations of $(p,k)$; Scenario a
is a worst case for the proposed SIS method.

| $(p,k)$ | MGP | CUSP | SIS
---|---|---|---|---
|  | Q0.5 | IQR | Q0.5 | IQR | Q0.5 | IQR
LPML |  00(16,4) |  0-28.68 | 0.42 |  0-28.68 | 0.43 |  0-28.65 | 0.41
|  00(32,8) |  0-60.08 | 0.45 |  0-60.09 | 0.45 |  0-60.07 | 0.49
|  0(64,12) | -117.68 | 0.56 | -117.75 | 0.53 | -117.88 | 0.56
| (128,16) | -225.04 | 1.04 | -225.13 | 1.04 | -228.76 | 1.47
$E(H_{a}\mid y)$ |  00(16,4) |  008.17 | 1.44 |  004.00 | 0.00 |  004.00 | 0.00
|  00(32,8) |  010.68 | 0.33 |  008.00 | 0.00 |  008.00 | 0.00
|  0(64,12) |  014.16 | 1.09 |  012.00 | 0.00 |  012.00 | 0.00
| (128,16) |  017.03 | 0.47 |  016.00 | 0.00 |  018.00 | 0.02

LPML, logarithm of the pseudo-marginal likelihood; CUSP, cumulative shrinkage
process; MGP, multiplicative gamma process; SIS, structured increasing
shrinkage process; Q0.5, median; IQR, interquartile range.

Scenario a is a worst case for the proposed method since there is no sparsity,
no structure, and the elements of the loadings matrix are similar in
magnitude. However, even in this case, structured increasing shrinkage
performs essentially identically to the best competitor, as illustrated by the
results in Table 1. The results of Ročková and George (2016) are not reported
as they are not competitive, as can be seen in table S2 in the Supplementary
Material. We report the median and interquartile range over the $R$ replicates
of the logarithm of the pseudo-marginal likelihood (Gelfand and Dey, 1994) and
of the estimated posterior mean of the number of factors $E(H_{a}\mid y)$.

Figure 2: Boxplots of mean squared error of the covariance matrix of each
model for different combinations of $(p,k,s)$ in Scenario b. Cov. MSE,
covariance mean squared error; CUSP, cumulative shrinkage process; MGP,
multiplicative gamma process; SIS, structured increasing shrinkage process.

Scenario b judges performance in detecting sparsity. The proposed approach
shows better performance in the logarithm of the pseudo-marginal likelihood
and mean squared error of the covariance matrix, particularly as sparsity
increases, as displayed in Fig. 2. Consistently with (Legramanti et al.,
2020), the covariance mean squared error is estimated in each simulation by
$\sum_{j,l}^{p}\sum_{t=1}^{S}(\omega_{jl}^{(t)}-\omega_{jl0})^{2}/\\{p(p+1)/2\\}$,
where $\omega_{jl0}$ and $\omega_{jl}^{(t)}$ are the elements $jl$ of
$\Omega_{0}=\Lambda_{0}\Lambda_{0}^{\mathrm{\scriptscriptstyle T}}+I_{p}$ and
$\Omega^{(t)}=\Lambda^{(t)}\Lambda^{(t){\mathrm{\scriptscriptstyle
T}}}+I_{p}$, respectively. The proposed approach allows exact zeros in the
loadings, while the competitors require thresholding to infer sparsity.
Following the thresholding approach described in Section S2.2 of the
Supplementary Material, we evaluate performance in inferring the sparsity
pattern via the mean classification error:

$MCE=\frac{1}{S}\sum_{t=1}^{S}\frac{\sum_{j=1}^{p}\sum_{h=1}^{k^{*(t)}}|\mathbbm{1}{(\lambda_{jh0}=0)}-\mathbbm{1}{(\lambda_{jh}^{(t)}=0)}|}{pk},$

where $k^{*(t)}$ is the maximum between the true number of factors $k$ and
$H_{a}^{(t)}$, and $\lambda_{jh0}$ and $\lambda_{jh}^{(t)}$ are the elements
$jh$ of $\Lambda_{0}$ and $\Lambda^{(t)}$, respectively. If $H_{a}^{(t)}$ or
$k$ are smaller than $k^{*}$, we fix the higher indexed columns at zero,
possibly leading to a mean classification error bigger than one. The results
reported in Table 2 show that the proposed structured increasing shrinkage
prior is much more effective in identifying sparsity in $\Lambda$, maintaining
good performance even with large $p$ and in strongly sparse contexts. Also,
more accurate estimation of the number of factors is obtained, as reported in
Table S1 in the Supplementary Material.

Table 2: Median and interquartile range of the mean classification error
computed in 25 replications assuming Scenario b and several combinations of
$(p,k,s)$

MCE | $(p,k,s)$ | MGP | CUSP | SIS
---|---|---|---|---
|  | Q0.5 | IQR | Q0.5 | IQR | Q0.5 | IQR
|  00(16,4,0.6) | 1.06 | 0.16 | 0.38 | 0.01 | 0.24 | 0.09
|  00(32,8,0.4) | 0.70 | 0.07 | 0.48 | 0.08 | 0.16 | 0.09
|  0(64,12,0.3) | 0.61 | 0.07 | 0.58 | 0.01 | 0.09 | 0.06
| (128,16,0.2) | 0.49 | 0.03 | 0.52 | 0.08 | 0.04 | 0.01

MCE, mean classification error; MGP, multiplicative gamma process; CUSP,
cumulative shrinkage process; SIS, structured increasing shrinkage process;
Q0.5, median; IQR, interquartile range.

Similar comments apply in Scenarios c and d reported in Fig. S2 in the
Supplementary Material. The superior performance of the structured increasing
shrinkage model is only partially mitigated in Scenario c for large $p$ for
the logarithm of the pseudo-marginal likelihood. In Scenario d, the use of
meta covariates has a mild benefit in identifying the sparsity pattern. In
lower signal-to-noise settings, meta covariates have a bigger impact, and they
also aid interpretation, as illustrated in the next section. Additional
details, tables, and plots for all scenarios are reported in Section S2.3 of
the Supplementary Material.

##  5 Finnish bird co-occurrence application

We illustrate our approach by modelling co-occurrence of the fifty most common
bird species in Finland (Lindström et al., 2015), focusing on data in 2014.
Response $y$ is an $n\times p$ binary matrix denoting occurrence of $p=50$
species in $n=137$ sampling areas. An $n\times c$ environmental covariate
matrix $w$ is available, including a 5-level habitat type, ‘spring
temperature’ (mean temperature in April and May), and $(\mbox{spring
temperature})^{2}$, leading to $c=7$. We consider a meta covariate $p\times q$
matrix $x$ of species traits: logarithm of typical body mass, migratory
strategy (short-distance migrant, resident species, long-distance migrant),
and a 7-level superfamily index. We model species presence or absence via a
multivariate probit regression model:

$y_{ij}=\mathbbm{1}{(z_{ij}>0)},\quad
z_{ij}=w_{i}^{T}\mu_{j}+\epsilon_{ij},\quad\epsilon_{i}=(\epsilon_{i1},\ldots,\epsilon_{ip})^{T}\sim
N_{p}(0,\Lambda\Lambda^{T}+I_{p}),$ (6)

where $\mu_{j}$ characterizes impact of environmental covariates on species
occurrence probabilities, and covariance in the latent $z_{i}$ vector is
characterized through a factor model. To borrow information across species
while incorporating species traits, we let

$\mu_{j}\sim N_{c}(b\,x_{j},\sigma^{2}_{\mu}I_{c}),\quad
b=(b_{1},\ldots,b_{q}),\quad b_{m}\sim N_{c}(0,\sigma^{2}_{b}I_{c}),$ (7)

where $b$ is a $c\times q$ coefficient matrix with column vectors $b_{m}$
given Gaussian priors.

Model (6)–(7) is consistent with popular joint species distribution models
(Ovaskainen et al., 2016; Tikhonov et al., 2017; Ovaskainen and Abrego, 2020),
with current standard practice using a multiplicative gamma process for
$\Lambda$. We compare this approach to an analysis that instead uses our
proposed structured increasing shrinkage prior to allow the species traits $x$
to impact $\Lambda$ and hence the covariance structure across species. After
standardizing $w$ and $x$, we set $\alpha=4$, $a_{\theta}=b_{\theta}=2$ and
$\sigma_{\mu}=\sigma_{b}=1$. Posterior sampling is straightforward via a Gibbs
sampler reported in Section S3.1 of the Supplementary Material.

Figure S8 in the Supplementary Material displays the posterior means of $\mu$
and $b$. A first investigation shows large heterogeneity of the habitat type
effects across different species. Matrix $b$ shows that covariate effects tend
to not depend on migratory strategy or body mass, with the exception of urban
habitats tending to have more migratory birds.

The estimated $\Lambda$ and meta covariate coefficients $\beta$, following the
guidelines of Section 3.3, are displayed in Fig. 3.

Figure 3: Posterior summaries $\Lambda^{(t^{*})}$ and $\beta^{(t^{*})}$ of the
structured increasing shrinkage model; rows of left matrix refer to 50 birds
species, and rows of right matrix to ten species traits. Light coloured cells
of $\beta^{(t^{*})}$ induce shrinkage on corresponding cells of
$\Lambda^{(t^{*})}$.

The loadings matrix is quite sparse, indicating that each latent factor
impacts a small group of species. Positive sign of the loadings means that
high levels of the corresponding factors increase the probability of observing
birds from those species. Lower elements of $\beta^{(t^{*})}$, represented
with light cells on the right panel, induce higher shrinkage on the
corresponding group of birds. To facilitate interpretation, we re-arrange the
rows of $\Lambda^{(t^{*})}$ according to the most relevant species traits in
terms of shrinkage, which are migration strategy and body mass. The species
influenced by the first factor are fairly homogeneous, characterized by short
distance or resident migratory strategies and a larger body mass. The strongly
negative value of $\beta_{(t^{*})42}$ suggests heavier species of birds tend
to have loadings close to zero for the second factor. This is also true for
the third factor, which also does not impact short-distance migrants.

Figure S9 in the Supplementary Material shows a spatial map of the sampling
units coloured accordingly to the values of the first and the third latent
factors. We can interpret these latent factors as unobserved environmental
covariates. We find that the species traits included in our analysis only
partially explain the loadings structure; this is as expected and provides
motivation for the proposed approach.

Figure 4: Posterior mean of the correlation matrices estimated by the
structured increasing shrinkage model (on the left) and the multiplicative
gamma process model (on the right).  Figure 5:  Graphical representation based
on the inverse of the posterior mean of the correlation matrices estimated by
the structured increasing shrinkage model (on the left) and the multiplicative
gamma process model (on the right). Edge thicknesses are proportional to the
latent partial correlations between species. Values below 0.025 are not
reported. Nodes are positioned using a Fruchterman–Reingold force-direct
algorithm.

Sparsity in the loadings matrix helps in interpretation. Species may load on
the same factor not just because they have similar traits but also because
they tend to favor similar habitats for reasons not captured by the measured
traits.

The induced covariance matrix $\Omega=\Lambda\Lambda^{T}+I_{p}$ across species
is of particular interest. We compare estimates of $\Omega$ under the
multiplicative gamma process, estimated using the R package hmsc (Tikhonov et
al., 2020), and our proposed structured increasing shrinkage model. Figure 4
reports the posterior mean of the correlation matrices under the two competing
models. The network graph based on the posterior mean of the partial
correlation matrices, reported in Fig. 5, reveals several communities of
species under the proposed structured increasing shrinkage prior that are not
evident under the multiplicative gamma.

We also find that the multiplicative gamma process provides a slightly worse
fit to the data. The logarithm of the pseudo marginal likelihood computed on
the posterior samples of the structured increasing shrinkage model is equal to
$-21.06$, higher than that achieved by the competing model, which is $-21.36$.
Using 4-fold cross-validation, we compared the log-likelihood evaluated in the
held-out data, with $\mu$ and $\Omega$ estimated by the posterior mean in the
training set. The mean of the log-likelihood was $-22.62$ under the structured
increasing shrinkage and $-23.22$ under the multiplicative gamma process
prior.

## Acknowledgement

The authors thank Daniele Durante, Sirio Legramanti, Otso Ovaskainen, and Gleb
Tikhonov for useful comments on an early version of this manuscript.

## Appendix

### Lemmas and proofs

###### of Theorem 1.

The variance of $\lambda_{jh}$ is

$\text{var}(\lambda_{jh})=E\\{E(\lambda_{jh}^{2}\mid\phi_{jh},\gamma_{h},\tau_{0})\\}=E\\{E(\theta_{jh}\mid\phi_{jh},\gamma_{h},\tau_{0})\\}.$

By construction,
$E(\theta_{jh}\mid\phi_{jh},\gamma_{h},\tau_{0})=\phi_{jh}\gamma_{h}\tau_{0}$.
Then,

$\text{var}(\lambda_{jh})=E(\phi_{jh}\gamma_{h}\tau_{0})=E(\phi_{j1})E(\gamma_{h})E(\tau_{0})>E(\phi_{j1})E(\gamma_{h+1})E(\tau_{0})=\text{var}(\lambda_{jh+1}),$

since the scale parameters are independent and the local scale $\phi_{jh}$ is
equally distributed over the column index $h$. ∎

To prove Proposition 2 we need to introduce the following Lemma.

######  Lemma 3.

Let $u,v$ denote two real positive random variables. If at least one among
$(u\mid v)$ and $(v\mid u)$ is power law tail distributed, then the product
$uv$ is power law tail distributed.

###### Proof.

For a positive value $w$, we can write

$\displaystyle\text{pr}(uv>w)$
$\displaystyle=\int_{0}^{\infty}\text{pr}(u>w/v\mid
v)f(v)\text{d}v=E\\{F_{u\mid v}^{C}(w/v)\\},$

where $F_{u\mid v}^{C}(w)=\text{pr}(u>w\mid v)$ and $f(v)$ is the probability
density function of $v$. If $F_{u\mid v}^{C}(w)\geq c\,w^{-\alpha}$ with
$c,\alpha$ positive constants and $w$ greater than a sufficiently large number
$L$, then

$\text{pr}(uv>w)\geq
E\\{c\,(w/v)^{-\alpha}\\}=c\,w^{-\alpha}E(v^{\alpha})\quad w>L\gg 0.$

If $E(v^{\alpha})=\infty$, then
$\text{pr}(uv>w)>c\,w^{-\alpha}=O(w^{-\alpha})$, otherwise
$\text{pr}(uv>w)\geq\nu(w)$ for $w>L$, with $\nu(w)$ a function of order
$O(w^{-\alpha})$ as $w$ goes to infinity. This shows that the right tail of
the distribution of the random variable $uv$ follows a power law behaviour. ∎

###### of Proposition 2.

Consider the strictly positive random variables
$\theta_{jh}^{*}=(\theta_{jh}\mid\theta_{jh}>0)$,
$\tau_{0}^{*}=(\tau_{0}\mid\tau_{0}>0)$,
$\gamma_{h}^{*}=(\gamma_{h}\mid\gamma_{h}>0)$, and
$\phi_{jh}^{*}=(\phi_{jh}\mid\phi_{jh}>0)$. Since $\theta_{jh}^{*}$ is equal
to the product $\tau^{*}\gamma_{h}^{*}\phi_{jh}^{*}$ of independent positive
random variables, Lemma 3 ensures that if at least one of those scale
parameters follows a power law tail distribution, then $\theta_{jh}^{*}$ is
power law tail distributed, so that $\text{pr}(\theta_{jh}^{*}>\theta)\geq
c\,\theta^{-\alpha}$ for $c,\alpha$ positive constants and $\theta>L$. Without
loss of generality, we focus on the right tail of $\lambda_{jh}$. Let

$\text{pr}(\lambda_{jh}>\lambda)=\text{pr}(\lambda_{jh}>\lambda\mid\theta_{jh}>0)\,\text{pr}(\theta_{jh}>0)+\text{pr}(\lambda_{jh}>\lambda\mid\theta_{jh}=0)\,\text{pr}(\theta_{jh}=0).$
(8)

It is straightforward to observe that $\lambda_{jh}$ marginally has a power
law tail if and only if $(\lambda_{jh}\mid\theta_{jh}>0)$ is power law tail
distributed and $\text{pr}(\theta_{jh}>0)$ is strictly positive. Since
$\text{pr}(\tau_{0}>0)>0$, $\text{pr}(\gamma_{h}>0)>0$, and
$\text{pr}(\phi_{jh}>0)>0$, then $\text{pr}(\theta_{jh}>0)>0,$ given
independence between the scale parameters. Focusing on $\theta_{jh}>0$ in the
first term of the right hand side of (8), we have

$\text{pr}(\lambda_{jh}>\lambda\mid\theta_{jh}^{*})=1-\Phi(\lambda\,\theta_{jh}^{*\,-0.5}),$

and we want to prove that the marginal
$F^{c}_{\lambda_{jh}}(\lambda)=\text{pr}(\lambda_{jh}>\lambda)$ is sub-
exponential as $\lambda\rightarrow\infty$. Using the lower bound for the right
tail of the standard Gaussian of Abramowitz and Stegun (1948),

$1-\Phi(\lambda\,\theta_{jh}^{*\,-0.5})\geq\left(\frac{2}{\pi}\right)^{0.5}\frac{\theta_{jh}^{*\,0.5}}{\lambda+(\lambda^{2}+4\theta_{jh}^{*})^{0.5}}e^{-\lambda^{2}/(2\theta_{jh}^{*})}.$

Marginalizing over $\theta_{jh}^{*}$, we obtain

$\text{pr}(\lambda_{jh}>\lambda\mid\theta_{jh}^{*})\geq
E\left\\{\left(\frac{2}{\pi}\right)^{0.5}\frac{\theta_{jh}^{*\,0.5}}{\lambda+(\lambda^{2}+4\theta_{jh}^{*})^{0.5}}e^{-\lambda^{2}/(2\theta_{jh}^{*})}\right\\}=E\left\\{t_{\lambda}(\theta_{jh}^{*})\right\\},$

where $t_{\lambda}(\theta_{jh}^{*})$ is a monotonically increasing nonnegative
function defined on the positive real line. Applying Markov’s inequality, we
have
$E\\{t_{\lambda}(\theta_{jh}^{*})\\}>\text{pr}(\theta_{jh}^{*}>\epsilon)t_{\lambda}(\epsilon)$,
and letting $\epsilon=\lambda^{2}$

$E\left\\{t_{\lambda}(\theta_{jh}^{*})\right\\}>\text{pr}(\theta_{jh}^{*}>\lambda^{2})\frac{e^{-0.5}}{1+5^{0.5}}\left(\frac{2}{\pi}\right)^{0.5}.$

If $\text{pr}(\theta_{jh}^{*}>\lambda)\geq c\,\lambda^{-\alpha}$ for certain
$\alpha,c$ positive constants and $\lambda$ sufficiently large, then

$\text{pr}(\lambda_{jh}>\lambda\mid\theta_{jh}^{*})\geq\frac{e^{-0.5}}{1+5^{0.5}}\left(\frac{2}{\pi}\right)^{0.5}c\,\lambda^{-2\alpha}=\tilde{c}\lambda^{-\tilde{\alpha}},$

where $\tilde{c}=e^{-0.5}(1+5^{0.5})^{-1}(2/\pi)^{0.5}c>0$ and
$\tilde{\alpha}=\alpha/2>0$. By symmetry,
$\text{pr}(\lambda_{jh}<-\lambda\mid\theta_{jh}>0)\geq\tilde{c}\lambda^{-\tilde{\alpha}}$
for $\lambda>L$ sufficiently large. It is sufficient that the marginal
distribution of $\theta_{jh}^{*}$ has power law right tail to guarantee that
$(\lambda_{jh}\mid\theta_{jh}>0)$ has power law tail and then that marginally
$\lambda_{jh}$ has power law tail. ∎

###### of Theorem 2.

The mode of the conditional posterior density of $\lambda_{jh}$ is
$\tilde{\lambda}_{jh}$ such that

$l_{s}(\tilde{\lambda}_{jh};y,\eta)+\frac{\partial}{\partial\lambda}\log\\{f_{\lambda_{jh}\mid\Lambda_{-jh}}(\lambda)\\}\big{|}_{\lambda=\tilde{\lambda}_{jh}}=0,$
(9)

where $l_{s}(\tilde{\lambda}_{jh};y,\eta)$ is the $jh$th element of the score
function of the likelihood for the data $y$ conditionally on the latent
variables $\eta$, and $f_{\lambda_{jh}\mid\Lambda_{-jh}}$ is the conditional
prior density function of $(\lambda_{jh}\mid\Lambda_{-jh})$. Given prior
symmetry, without loss of generality, we focus on $\hat{\lambda}_{jh}>0$. In a
neighbourhood
$(\hat{\lambda}_{jh}-\varepsilon,\hat{\lambda}_{jh}+\varepsilon)$ of the
conditional maximum likelihood estimate $\hat{\lambda}_{jh}$ of
$\lambda_{jh}$, we can approximate the score function using a Taylor
expansion:
$l_{s}(\lambda;y)=-\mathcal{J}(\hat{\lambda}_{jh})\,(\lambda-\hat{\lambda}_{jh})+\epsilon_{\varepsilon}$,
where $\mathcal{J}(\hat{\lambda}_{jh})>0$ is the negative of the derivative of
$l_{s}(\lambda;y)$ evaluated at $\lambda=\hat{\lambda}_{jh}$, and
$\epsilon_{\varepsilon}$ is an approximation error term such that
$\lim_{\varepsilon\to 0}\epsilon_{\varepsilon}/\varepsilon=0$. For
$\hat{\lambda}_{jh}$ large enough, such that
$\hat{\lambda}_{jh}-\varepsilon>L$ with $L\gg 0$, Lemma 1 holds for every
$\lambda$ in
$(\hat{\lambda}_{jh}-\varepsilon,\hat{\lambda}_{jh}+\varepsilon)$, leading to
the lower bound

$-\mathcal{J}(\hat{\lambda}_{jh})\,(\lambda-\hat{\lambda}_{jh})+f_{lb}^{\prime}(\lambda)+\epsilon_{\varepsilon}\leq
l_{s}(\lambda;y)+\frac{\partial}{\partial\lambda}\log\\{f_{\lambda_{jh}\mid\Lambda_{-jh}}(\lambda)\\},$

where $f_{lb}^{\prime}(\lambda)$ is a non positive continuous function for
every $\lambda>0$, $\lim_{\lambda\to+\infty}f_{lb}^{\prime}(\lambda)=0$. Let
$\varepsilon$ be a function of $\hat{\lambda}_{jh}$ such that
$\lim_{\hat{\lambda}_{jh}\to\infty}\varepsilon=0$ and
$\lim_{\hat{\lambda}_{jh}\to\infty}f_{lb}^{\prime}(\hat{\lambda}_{jh})/\varepsilon=0.$
The limit for $\hat{\lambda}_{jh}\to\infty$ of the lower bound evaluated in
$\hat{\lambda}_{jh}-\varepsilon$ is

$\lim_{\hat{\lambda}_{jh}\to\infty}\mathcal{J}(\hat{\lambda}_{jh})\,\varepsilon+f_{lb}^{\prime}(\hat{\lambda}_{jh}-\varepsilon)+\epsilon_{\varepsilon}=\lim_{\hat{\lambda}_{jh}\to\infty}|\varepsilon|\,\\{\mathcal{J}(\hat{\lambda}_{jh})+f_{lb}^{\prime}(\hat{\lambda}_{jh}-\varepsilon)/|\varepsilon|+\epsilon_{\varepsilon}/|\varepsilon|\\}.$

Under Assumption 1,
$\lim_{\hat{\lambda}_{jh}\to\infty}\mathcal{J}(\hat{\lambda}_{jh})+f_{lb}^{\prime}(\hat{\lambda}_{jh}-\varepsilon)/|\varepsilon|+\epsilon_{\varepsilon}/|\varepsilon|\geq
0,$ which guarantees
$\hat{\lambda}_{jh}-\varepsilon\leq\tilde{\lambda}_{jh}\leq\hat{\lambda}_{jh},$
and, hence
$\lim_{\hat{\lambda}_{jh}\to\infty}|\tilde{\lambda}_{jh}-\hat{\lambda}_{jh}|=0,$
which proves the theorem. ∎

###### of Theorem 3.

Since the local scales are independent, conditionally on $\beta$, we can apply
the Chernoff’s method and obtain the following upper bound

$\text{pr}\\{|\text{supp}_{\epsilon}(\lambda_{h})|>as_{p}\mid\beta_{h},\gamma_{h},\tau_{0}\\}\leq\exp(-s_{p}a\,t)\,\exp\bigg{\\{}(e^{t}-1)\,\sum_{j=1}^{p}\zeta_{\epsilon
jh}\bigg{\\}},$

for every $t>0$ and $\zeta_{\epsilon
jh}=\\{\tau_{0}\,\gamma_{h}\,g(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})\\}/\epsilon^{2}$ a function of $\beta_{h}$. Since
$g(x_{j}^{\mathrm{\scriptscriptstyle T}}\beta_{h})$ is of order $\leq
O(\log(p)/p)$ by assumption and is limited above with respect to $\beta_{h}$,
we can deduce $g(x_{j}^{\mathrm{\scriptscriptstyle T}}\beta_{h})\leq
c_{j}\log(p)/p$ for $p$ sufficiently large and for some constant $c_{j}>0$
that does not depend on $\beta_{h}$ and is asymptotically of order $O(1)$ with
respect to $p$. Then, for $p\gg 0$,

$\displaystyle\sum_{j=1}^{p}g(x_{j}^{\mathrm{\scriptscriptstyle T}}\beta_{h})$
$\displaystyle\leq\sum_{j=1}^{p}c_{j}\log(p)/p\leq p\,\log(p)/p\,\max_{1\leq
j\leq p}c_{j}=m\log(p),$

where $m=\max_{1\leq j\leq p}c_{j}$ does not depend on $\beta_{h}$. Then, the
upper bound is

$\text{pr}\\{|\text{supp}_{\epsilon}(\lambda_{h})|>as_{p}\mid\beta_{h},\gamma_{h},\tau_{0}\\}\leq\exp\bigg{\\{}-s_{p}a\,t+(e^{t}-1)\,\frac{\tau_{0}\gamma_{h}}{\epsilon^{2}}\,m\log(p)\bigg{\\}}.$

Let us choose $t=\log\\{\epsilon^{2}/(\tau_{0}\gamma_{h}m)+1\\}$. Since
$s_{p}\geq\log(p)\,c_{s}$ for a certain $c_{s}>0$, then, for any
$a>(c_{s}t)^{-1}$, we can write

$\text{pr}\\{|\text{supp}_{\epsilon}(\lambda_{h})|>as_{p}\mid\beta_{h},\gamma_{h},\tau_{0}\\}\leq\exp\bigg{\\{}-\log(p)\,\tilde{a}\bigg{\\}},$

where $\tilde{a}$ is a positive constant such that
$a=(1+\tilde{a})(c_{s}t)^{-1}$. The upper bound does not depend on
$\beta_{h}$, so

$\text{pr}\\{|\text{supp}_{\epsilon}(\lambda_{h})|>as_{p}\mid\gamma_{h},\tau_{0}\\}\leq\nu(p)$

with $\nu(p)$ of order $O(p^{-1})$ that goes to zero. ∎

## References

  * Abramowitz and Stegun (1948) Abramowitz, M. and Stegun, I. A. (1948).  Handbook of mathematical functions with formulas, graphs, and mathematical tables, volume 55.  US Government printing office.
  * An et al. (2013) An, X., Yang, Q., and Bentler, P. M. (2013).  A latent factor linear mixed model for high-dimensional longitudinal data analysis.  Statistics in medicine, 32(24):4229–4239.
  * Aßmann et al. (2016) Aßmann, C., Boysen-Hogrefe, J., and Pape, M. (2016).  Bayesian analysis of static and dynamic factor models: An ex-post approach towards the rotation problem.  Journal of Econometrics, 192(1):190–206.
  * Bhattacharya and Dunson (2011) Bhattacharya, A. and Dunson, D. B. (2011).  Sparse Bayesian infinite factor models.  Biometrika, 98(2):291–306.
  * Bhattacharya et al. (2015) Bhattacharya, A., Pati, D., Pillai, N. S., and Dunson, D. B. (2015).  Dirichlet-Laplace priors for optimal shrinkage.  Journal of the American Statistical Association, 110(512):1479–1490.
  * Carvalho et al. (2010) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2010).  The horseshoe estimator for sparse signals.  Biometrika, 97(2):465–480.
  * Dahl (2006) Dahl, D. B. (2006).  Model-based clustering for expression data via a Dirichlet process mixture model.  Bayesian inference for gene expression and proteomics, 4:201–218.
  * Durante (2017) Durante, D. (2017).  A note on the multiplicative gamma process.  Statistics and Probability Letters, 122:198–204.
  * Ferrari and Dunson (2020) Ferrari, F. and Dunson, D. B. (2020).  Bayesian factor analysis for inference on interactions.  Journal of the American Statistical Association, pages 1–12.
  * Gelfand and Dey (1994) Gelfand, A. E. and Dey, D. K. (1994).  Bayesian model choice: asymptotics and exact calculations.  Journal of the Royal Statistical Society: Series B (Methodological), 56(3):501–514.
  * Jun and Tao (2013) Jun, L. and Tao, D. (2013).  Exponential Family Factors for Bayesian Factor Analysis.  IEEE Transactions on neural networks and learning systems, 24(6):964—-976.
  * Legramanti et al. (2020) Legramanti, S., Durante, D., and Dunson, D. B. (2020).  Bayesian cumulative shrinkage for infinite factorizations.  Biometrika, 107(3):745–752.
  * Lindström et al. (2015) Lindström, Å., Green, M., Husby, M., Kålås, J. A., and Lehikoinen, A. (2015).  Large-scale monitoring of waders on their boreal and arctic breeding grounds in northern Europe.  Ardea, 103(1):3–15.
  * Liu and Vandenberghe (2010) Liu, Z. and Vandenberghe, L. (2010).  Interior-point method for nuclear norm approximation with application to system identification.  SIAM Journal on Matrix Analysis and Applications, 31(3):1235–1256.
  * Lopes and West (2004) Lopes, H. F. and West, M. (2004).  Bayesian model assessment in factor analysis.  Statistica Sinica, pages 41–67.
  * McParland et al. (2014) McParland, D., Gormley, I. C., McCormick, T. H., Clark, S. J., Kabudula, C. W., and Collinson, M. A. (2014).  Clustering south African households based on their asset status using latent variable models.  The annals of applied statistics, 8(2):747.
  * Miller et al. (2019) Miller, J. E., Li, D., LaForgia, M., and Harrison, S. (2019).  Functional diversity is a passenger but not driver of drought-related plant diversity losses in annual grasslands.  Journal of Ecology, 107(5):2033–2039.
  * Miller and Harrison (2018) Miller, J. W. and Harrison, M. T. (2018).  Mixture models with a prior on the number of components.  Journal of the American Statistical Association, 113(521):340–356.
  * Mitchell and Beauchamp (1988) Mitchell, T. J. and Beauchamp, J. J. (1988).  Bayesian variable selection in linear regression.  Journal of the American Statistical Association, 83(404):1023–1036.
  * Mnih and Salakhutdinov (2008) Mnih, A. and Salakhutdinov, R. R. (2008).  Probabilistic matrix factorization.  In Advances in neural information processing systems, pages 1257–1264.
  * Montagna et al. (2012) Montagna, S., Tokdar, S. T., Neelon, B., and Dunson, D. B. (2012).  Bayesian latent factor regression for functional and longitudinal data.  Biometrics, 68(4):1064–1073.
  * Murray et al. (2013) Murray, J. S., Dunson, D. B., Carin, L., and Lucas, J. E. (2013).  Bayesian Gaussian copula factor models for mixed data.  Journal of the American Statistical Association, 108(502):656–665.
  * Ovaskainen and Abrego (2020) Ovaskainen, O. and Abrego, N. (2020).  Joint Species Distribution Modelling: With Applications in R.  Cambridge University Press.
  * Ovaskainen et al. (2016) Ovaskainen, O., Abrego, N., Halme, P., and Dunson, D. (2016).  Using latent variable models to identify large networks of species-to-species associations at different spatial scales.  Methods in Ecology and Evolution, 7(5):549–555.
  * Polson and Scott (2010) Polson, N. G. and Scott, J. G. (2010).  Shrink globally, act locally: Bayesian sparsity and regularization.  Bayesian Statistics, 9:1–16.
  * Polson et al. (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using pólya–gamma latent variables.  Journal of the American statistical Association, 108(504):1339–1349.
  * Reich and Bandyopadhyay (2010) Reich, B. J. and Bandyopadhyay, D. (2010).  A latent factor model for spatial data with informative missingness.  The annals of applied statistics, 4(1):439.
  * Roberts and Rosenthal (2007) Roberts, G. O. and Rosenthal, J. S. (2007).  Coupling and ergodicity of adaptive Markov chain Monte Carlo algorithms.  Journal of applied probability, 44(2):458–475.
  * Ročková and George (2016) Ročková, V. and George, E. I. (2016).  Fast bayesian factor analysis via automatic rotations to sparsity.  Journal of the American Statistical Association, 111(516):1608–1622.
  * Rousseau and Mengersen (2011) Rousseau, J. and Mengersen, K. (2011).  Asymptotic behaviour of the posterior distribution in overfitted mixture models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(5):689–710.
  * Roweis and Ghahramani (1999) Roweis, S. and Ghahramani, Z. (1999).  A unifying review of linear Gaussian models.  Neural computation, 11(2):305–345.
  * Roy et al. (2019) Roy, A., Schaich-Borg, J., and Dunson, D. B. (2019).  Bayesian time-aligned factor analysis of paired multivariate time series.  arXiv preprint arXiv:1904.12103.
  * Schiavon and Canale (2020) Schiavon, L. and Canale, A. (2020).  On the truncation criteria in infinite factor models.  Stat, 9:e298.
  * Thomas et al. (2009) Thomas, D. C., Conti, D. V., Baurley, J., Nijhout, F., Reed, M., and Ulrich, C. M. (2009).  Use of pathway information in molecular epidemiology.  Human genomics, 4(1):21.
  * Tikhonov et al. (2017) Tikhonov, G., Abrego, N., Dunson, D., and Ovaskainen, O. (2017).  Using joint species distribution models for evaluating how species-to-species associations depend on the environmental context.  Methods in Ecology and Evolution, 8(4):443–452.
  * Tikhonov et al. (2020) Tikhonov, G., Opedal, Ø. H., Abrego, N., Lehikoinen, A., de Jonge, M. M., Oksanen, J., and Ovaskainen, O. (2020).  Joint species distribution modelling with the R-package hmsc.  Methods in ecology and evolution, 11(3):442–447.
  * Wade et al. (2018) Wade, S., Ghahramani, Z., et al. (2018).  Bayesian cluster analysis: Point estimation and credible balls (with discussion).  Bayesian Analysis, 13(2):559–626.
  * Yang et al. (2018) Yang, L., Fang, J., Duan, H., Li, H., and Zeng, B. (2018).  Fast low-rank Bayesian matrix completion with hierarchical Gaussian prior models.  IEEE Transactions on Signal Processing, 66(11):2804–2817.
  * Yuan and Lin (2006) Yuan, M. and Lin, Y. (2006).  Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(1):49–67.

Supplementary Material for "Generalized infinite factorization models"

This supplementary material available at Biometrika online includes the
statement and proof of Proposition S1 and the proofs of Proposition 1, Lemmas
1–2, and Corollaries 1–3. The Gibbs sampling algorithm, settings, and
additional results of the simulations and ecology data analysis are reported,
including trace plots and a sensitivity analysis to varying hyperparameters.

##  Appendix S1 Propositions and proofs

######  Proposition S3.

Let $\Pi_{\Lambda}\otimes\Pi_{\Sigma}$ denote the prior on $(\Lambda,\Sigma)$.
Let $\Theta_{\Lambda}$ and $\Theta_{\Sigma}$ denote the sample spaces of the
matrices $\Lambda$ and $\Sigma$, respectively. If
$\sum_{h=1}^{\infty}E(\gamma_{h})<\infty$, then,
$\Pi_{\Lambda}\otimes\Pi_{\Sigma}(\Theta_{\Lambda}\times\Theta_{\Sigma})=1$.

###### Proof.

Assume $\Sigma\in\Theta_{\Sigma}$ and
$(\Psi,\Lambda)\in\Theta_{\Psi}\times\Theta_{\Lambda}$, with $\Theta_{\Sigma}$
the set of $p\times p$ positive semi-definite matrices with finite elements,
and

$\displaystyle\Theta_{\Psi}\times\Theta_{\Lambda}$
$\displaystyle=\bigg{\\{}\Lambda=(\lambda_{jh}),\,\Psi=(\psi_{hh}):\sum_{h=1}^{\infty}\lambda_{jh}\psi_{hh}\lambda_{sh}<\infty\;\forall\
j,s\in(1,\ldots,p)\bigg{\\}}.$

Due to independence, we can study the prior on $\Sigma$ and $\Lambda$
separately. The prior on $\Sigma$ is defined on the set of positive semi-
definite matrices. Therefore, it is sufficient to prove that the elements of
$\Lambda\Psi\Lambda^{\mathrm{\scriptscriptstyle T}}$ are finite almost surely.
Using Cauchy-Schwartz, it is straightforward to show that all the entries of
$\Lambda\Lambda^{\mathrm{\scriptscriptstyle T}}$ are finite if and only if
$\sum_{h=1}^{\infty}\psi_{hh}\lambda_{jh}^{2}<\infty\ (j=1,\ldots,p)$. Let $c$
satisfy $c>\max_{h=1,\ldots,\infty}\psi_{hh}$ and let
$X_{j}^{{\mathrm{\scriptscriptstyle T}}}$ denote the $j$-th row vector of
matrix $X$. Since

$\displaystyle E(\lambda_{jh}^{2})$
$\displaystyle=E\\{E(\lambda_{jh}^{2}\mid\phi_{jh},\gamma_{h},\tau_{0})\\}=E(\phi_{jh})E(\gamma_{h})E(\tau_{0}),$
$\displaystyle E(\phi_{jh})$
$\displaystyle=E\\{E(\phi_{jh}\mid\beta_{h})\\}=E\\{g(X_{j}^{{\mathrm{\scriptscriptstyle
T}}}\beta_{h})\\}=E(\phi_{j1}),$

it is sufficient that $\sum_{h=1}^{\infty}E(\gamma_{h})<\infty$ to prove that
$\sum_{h=1}^{\infty}E(\lambda_{jh}^{2})=E(\phi_{j1})E(\tau_{0})\sum_{h=1}^{\infty}E(\gamma_{h})<\infty$
and then
$\sum_{h=1}^{\infty}\psi_{hh}\lambda_{jh}^{2}<c\sum_{h=1}^{\infty}\lambda_{jh}^{2}<\infty$.
∎

###### of Proposition 1.

The trace of $\Omega$ is
$\text{tr}({\Sigma})+\text{tr}(\Lambda_{H}\Psi_{H}\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}})+\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})$, where $\Lambda_{\Delta_{H}}=\Lambda-\Lambda_{H}$ and
$\Psi_{\Delta_{H}}=\Psi-\Psi_{H}$. Hence, it is equivalent to rewrite the
probability of interest as

$\text{pr}\left\\{\frac{\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})}{\text{tr}({\Omega})}\geq 1-T\right\\}.$

By Markov’s Inequality

$\text{pr}\left\\{\frac{\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})}{\text{tr}({\Omega})}\geq 1-T\right\\}\leq
E\left\\{\frac{\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})}{\text{tr}({\Omega})}\right\\}\big{/}(1-T).$

The expected ratio of two random variables $u$ and $v$ is
$E(u/v)=\text{cov}(u,1/v)+E(u)E(1/v)$, which allows us to write $E(u/v)\leq
E(u)E(1/v)$ if $\text{cov}(u,1/v)\leq 0$. Then, since the covariance between
$\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})$ and $\text{tr}({\Omega})$ is non-negative, the following inequality
holds

$E\left\\{\frac{\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})}{\text{tr}({\Omega})}\right\\}\leq
E\\{\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})\\}E\left(\frac{1}{\text{tr}({\Omega})}\right).$

The trace
$\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})$ is equal to
$\sum_{j=1}^{p}\sum_{h=H+1}^{\infty}\psi_{hh}\lambda_{jh}^{2}$. The variance
of $\lambda_{jh}$ is
$E(\lambda_{jh}^{2})=E(\phi_{j1})E(\gamma_{h})E(\tau_{0})$. Let $c$ satisfy
$c\geq\max_{h=1,\ldots,\infty}\psi_{hh}$. Since $E(\phi_{j1})$ is finite and
$E(\gamma_{h})=ab^{h-1}$ with $a,b$ positive constants and $b<1$, then

$E\\{\text{tr}(\Lambda_{\Delta_{H}}\Psi_{\Delta_{H}}\Lambda_{\Delta_{H}}^{\mathrm{\scriptscriptstyle
T}})\\}\leq cE(\tau_{0})a\frac{b^{H}}{1-b}\sum_{j=1}^{p}E(\phi_{j1}).$

Since
$\text{tr}({\Omega})=\text{tr}(\Lambda\Psi\Lambda^{\mathrm{\scriptscriptstyle
T}})+\text{tr}(\Sigma)$, we know that
$\text{tr}({\Omega})\geq\sum_{h=1}^{\infty}\psi_{hh}\lambda_{jh}^{2}+\sigma_{j}^{2}$
for any $j$ in $1,\ldots,p$, where $\sigma_{j}^{2}$ is the $j$-th diagonal
element of $\Sigma$. Then, for any $j$ in $1,\ldots,p$, we obtain

$\frac{1}{\text{tr}({\Omega})}\leq\frac{1}{\sum_{h=1}^{\infty}\psi_{hh}\lambda_{jh}^{2}+\sigma_{j}^{2}},$

and, consequently,

$E\left\\{\frac{1}{\text{tr}({\Omega})}\right\\}\leq
E\left(\sigma_{j}^{-2}\right),\qquad
E\left\\{\frac{1}{\text{tr}({\Omega})}\right\\}\leq
E\left\\{\left(\sum_{h=1}^{\infty}\psi_{hh}\lambda_{jh}^{2}\right)^{-1}\right\\}.$

Therefore, since
$m_{\Omega}=\min_{j=1,\ldots,p}\left[E(\sigma_{j}^{-2}),E\left\\{\left(\sum_{h=1}^{\infty}\psi_{hh}\lambda_{jh}^{2}\right)^{-1}\right\\}\right]<\infty$,
then

$\text{pr}\left\\{\frac{\text{tr}(\Lambda_{H}\Psi_{H}\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}})+\text{tr}(\Sigma)}{\text{tr}({\Omega})}\leq
T\right\\}\leq\bigg{(}\frac{1}{1-T}\bigg{)}a\,c\,\frac{b^{H}}{1-b}\,m_{\Omega}\,E(\tau_{0})\,\sum_{j=1}^{p}E(\phi_{j1}),$

as stated by the Theorem. ∎

###### of Lemma 1.

Consistently with Proposition 2, $(\lambda_{jh}\mid\Lambda_{-jh})$ has power
law tail if $(\theta_{jh}\mid\Lambda_{-jh})$ has power law tail. Furthermore,
$\text{pr}(|\lambda_{jh}|>\lambda\mid\Lambda_{-jh})$ has power law tail for
large $\lambda$ if and only if
$\text{pr}(|\lambda_{jh}|>\lambda\mid\Lambda_{-jh},\theta_{jh}>0)$ has power
law tail and $\text{pr}(\theta_{jh}>0\mid\Lambda_{-jh})>0$. The latter
inequality is always true when the marginal $\text{pr}(\theta_{jh}>0)$ is
positive.To prove $(\theta_{jh}\mid\Lambda_{-jh})$ has power law tail, we
apply Lemma 3. We first focus on proving the lemma when $\phi_{jh}$ satisfies
the power law tail condition with $\tau_{0}\gamma_{h}=w_{h}.$ As the local
scale $\phi_{jh}$ is independent from $(\Lambda_{-jh},w_{h})$ given
$\beta_{h}$, its conditional density is

$f_{\phi_{jh}\mid
w_{h},\Lambda_{-(jh)}}(\phi)=\int_{\Re}f_{\phi_{jh}\mid\beta_{h}}(\phi)\,f_{\beta_{h}\mid
w_{h},\Lambda_{-jh}}(\beta)\,\text{d}\beta.$

As the tail conditions hold for any possible prior on $\beta$, we have

$f_{\phi_{jh}}(\phi)=\int_{\Re}f_{\phi_{jh}\mid\beta_{h}}(x)\,f(\beta)\,\text{d}\beta,\qquad
f_{\phi_{jh}}(\tilde{\phi})\propto\tilde{\phi}^{-\alpha},\qquad\tilde{\phi}=\\{\phi:\phi>l\\},\qquad
L\gg 0,$

for any prior density $f$ defined on $\Re$. Hence, $(\phi_{jh}\mid
w_{jh},\Lambda_{-jh})$ is power law tail distributed. We now focus on proving
the lemma when $\tau_{0}$ or $\gamma_{h}$ are power law tail distributed. Let
$r_{h}^{*}=(r_{h}\mid r_{h}>0)$ and $w_{jh}^{*}=(w_{jh}\mid w_{jh}>0)$, where
$r_{h}$ is the scale parameter with power law tail and $w_{h}$ is the product
of the remaining two scale parameters, respectively. By Bayes’ Theorem

$f_{r_{h}^{*}\mid w_{jh}^{*},\Lambda_{-jh}}(r)=\frac{f_{\Lambda_{-jh}\mid
w_{jh}^{*},r_{h}^{*}}(\Lambda_{-jh};r)f_{r_{h}^{*}\mid
w_{jh}^{*}}(r)}{f_{\Lambda_{-jh}\mid w_{jh}^{*}}(\Lambda_{-jh})}.$

Since $r_{h}^{*}$ is independent from $w_{jh}^{*}$ for any parameter scale, it
is sufficient to prove that the function $f_{\Lambda_{-jh}\mid
w_{jh}^{*},r_{h}^{*}}(\Lambda_{-jh};r)$ decreases slower than
$c\,r^{\,-\alpha}$, for $c,\alpha>$ positive constants, when
$r\rightarrow\infty$. Denoting
$F_{\tau_{0},\phi_{11}\ldots\phi_{pk},\gamma_{1},\ldots,\gamma_{h-1},\gamma_{h+1},\ldots,\gamma_{k}\mid
w_{jh}^{*},r_{h}^{*}}$ the probability measure for conditional density
$f_{\tau_{0},\phi_{11}\ldots\phi_{pk},\gamma_{1},\ldots,\gamma_{h-1},\gamma_{h+1},\ldots,\gamma_{k}\mid
w_{jh}^{*},r_{h}^{*}}$, we can write

$\displaystyle f_{\Lambda_{-jh}\mid w_{jh}^{*},r_{h}^{*}}(\Lambda_{-jh};r)$
$\displaystyle=\int
f_{\Lambda_{-jh}\mid\tau_{0},\phi_{11}\ldots\phi_{pk},\gamma_{1},\ldots,\gamma_{k}}(\Lambda_{-jh};r)\,\text{d}F_{\tau_{0},\phi_{11}\ldots\phi_{pk},\gamma_{1},\ldots,\gamma_{h-1},\gamma_{h+1},\ldots,\gamma_{k}\mid
w_{jh}^{*},r_{h}^{*}}$
$\displaystyle=\int\prod_{(s,m)\neq(j,h)}f_{\lambda_{sm}\mid\theta_{sm}}(\lambda_{sm};r)\,\text{d}F_{\tau_{0},\phi_{11},\ldots,\gamma_{k}\mid
w_{jh}^{*},r_{h}^{*}}$
$\displaystyle=E\bigg{\\{}\prod_{(s,m)\neq(j,h)}f_{\lambda_{sm}\mid\theta_{sm}}(\lambda_{sm};r)\biggm{|}w_{jh}^{*},r_{h}^{*},\Lambda_{-jh}\bigg{\\}}$

The product inside the expectation is zero when there is a pair of indices
$(s,m)$ such that $\lambda_{sm}\neq 0$ and $\theta_{sm}=0$. However, since the
probability $\text{pr}(\theta_{sm}=0\mid\lambda_{sm}\neq 0)=0$, we know that
the expected value of the product between the functions
$f_{\lambda_{sm}\mid\theta_{sm}}(\lambda_{sm};r)$, given
$w_{jh}^{*},r_{h}^{*},\Lambda_{-jh}$, is strictly positive. We first focus on
the case $\gamma_{h}=r_{h}$ and prove that $f_{\Lambda_{-jh}\mid
w_{jh}^{*},\gamma_{h}^{*}}(\Lambda_{-jh};\gamma)$ decreases slower than
$c\gamma^{-\alpha}$ for $c,\alpha>0$. In this case, we can write the above
expectation as

$E\bigg{\\{}\prod_{s=1,m\neq h}^{p}f_{k}(\lambda_{sm})\,\prod_{s\neq
j}f_{\lambda_{sh}\mid\theta_{sh}}(\lambda_{sh};\gamma_{h}^{*})\biggm{|}w_{jh}^{*},\gamma_{h}^{*},\Lambda_{-jh}\bigg{\\}},$

where $\prod_{s=1,m\neq h}^{p}f_{k}(\lambda_{sm})$ is a product between
$(k-1)\times p$ strictly positive random variables that does not depend on
$w_{jh}^{*}$ and $\gamma_{h}^{*}$, while $\prod_{s\neq
j}f_{\lambda_{sh}\mid\theta_{sh}}(\lambda_{sh};\gamma_{h}^{*})$ is a product
between $p$ strictly positive random variables. In particular, if $w_{sh}=0$,
then
$f_{\lambda_{sm}\mid\theta_{sh}}(\lambda_{sh};\gamma_{h}^{*})=\mathbbm{1}{(\lambda_{sh}=0)}$.
If $w_{sh}>0$, then

$f_{\lambda_{sh}\mid\theta_{sh}}(\lambda_{sh};\gamma_{h}^{*})=(2\pi
w_{sh}^{*}\gamma_{h}^{*})^{-0.5}\exp\bigg{(}-\frac{\lambda_{sh}^{2}}{2w_{sh}^{*}\gamma_{h}^{*}}\bigg{)}>(2\pi
w_{sh}^{*}\gamma_{h}^{*})^{-0.5}\exp\bigg{(}-\frac{\lambda_{sh}^{2}}{2w_{sh}^{*}}\bigg{)}.$

Therefore, the upper bound

$f_{\lambda_{sh}\mid\theta_{sh}}(\lambda_{sh};\gamma_{h}^{*})\geq\begin{cases}\min\\{1,\,(2\pi
w_{sh}^{*}\gamma_{h}^{*})^{-0.5}\\},\qquad&\text{if}\;\lambda_{sh=0}\\\ (2\pi
w_{sh}^{*}\gamma_{h}^{*})^{-0.5}\exp\\{-\lambda_{sh}^{2}/(2w_{sh}^{*})\\}\qquad&\text{if}\;\lambda_{sh\neq
0},\end{cases}$

holds with probability equal to 1. For $\gamma>1$, we note that
$f_{\lambda_{sh}\mid\theta_{sh}}(\lambda_{sh};\gamma)\geq\gamma^{-0.5}u_{\lambda_{sh}}$
with

$u_{\lambda_{sh}}=\begin{cases}\min\\{1,\,(2\pi
w_{sh}^{*})^{-0.5}\\},\qquad&\text{if}\;\lambda_{sh}=0,\\\ (2\pi
w_{sh}^{*})^{-0.5}\exp\\{-\lambda_{sh}^{2}/(2w_{sh}^{*})\\}\qquad&\text{if}\;\lambda_{sh}\neq
0.\end{cases}$

Then,

$\displaystyle
E\bigg{\\{}\prod_{(s,m)\neq(j,h)}f_{\lambda_{sm}\mid\theta_{sm}}(\lambda_{sm};\gamma_{h}^{*})$
$\displaystyle\biggm{|}w_{jh}^{*},\gamma_{h}^{*},\Lambda_{-jh}\bigg{\\}}\geq$
$\displaystyle E\bigg{\\{}\prod_{s=1,m\neq
h}^{p}f_{k}(\lambda_{sm})\,\prod_{s\neq
j}\gamma_{h}^{*-0.5}u_{\lambda_{sh}}\biggm{|}w_{jh}^{*},\gamma_{h}^{*},\Lambda_{-jh}\bigg{\\}}=$
$\displaystyle\gamma_{h}^{*-0.5\,(p-1)}E\bigg{\\{}\prod_{s=1,m\neq
h}^{p}f_{k}(\lambda_{sm})\,\prod_{s\neq
j}u_{\lambda_{sh}}\biggm{|}w_{jh}^{*},\Lambda_{-jh}\bigg{\\}},$

where the expectation is strictly positive and not depending on $\gamma_{h}$.
Therefore, for $\gamma$ sufficiently large, $f_{\Lambda_{-jh}\mid
w_{jh}^{*},\gamma_{h}^{*}}(\Lambda_{-jh};\gamma)\geq c\gamma^{-\alpha},\qquad
c,\alpha>0$ holds, so that $(\gamma_{h}\mid w_{jh},\Lambda_{-jh})$ is power
law tail distributed. Similarly, if $\tau_{0}=r_{h}$ ($h=1,\ldots,H$),

$\displaystyle
E\bigg{\\{}\prod_{(s,m)\neq(j,h)}f_{\lambda_{sm}\mid\theta_{sm}}(\lambda_{sm};\tau_{0}^{*})$
$\displaystyle\mid w_{jh}^{*},\tau_{0}^{*},\Lambda_{-jh}\bigg{\\}}\geq$
$\displaystyle
E\bigg{(}\prod_{(s,m)\neq(j,h)}\tau_{0}^{*-0.5}u_{\lambda_{sm}}\mid
w_{jh}^{*},\tau_{0}^{*},\Lambda_{-jh}\bigg{)}=$
$\displaystyle\tau_{0}^{*-0.5\,(pH-1)}E\bigg{(}\prod_{(s,m)\neq(j,h)}u_{\lambda_{sm}}\mid
w_{jh}^{*},\Lambda_{-jh}\bigg{)},$

where

$u_{\lambda_{sm}}=\begin{cases}\min\\{1,\,(2\pi
w_{sm}^{*})^{-0.5}\\},\qquad\qquad&\text{if}\;\lambda_{sm}=0,\\\ (2\pi
w_{sm}^{*})^{-0.5}\exp\\{-\lambda_{sm}^{2}/(2w_{sm}^{*})\\},\qquad&\text{if}\;\lambda_{sm}\neq
0.\end{cases}$

is strictly positive and does not depend on $\tau_{0}$. Then, if the number
$H$ of columns of $\Lambda_{H}$ is finite, $f_{\Lambda_{-jh}\mid
w_{jh}^{*},\tau_{0}^{*}}(\Lambda_{-jh};\tau)\geq c\tau^{-\alpha}$ with
$c,\alpha>0$ and $\tau$ sufficiently large, implying $(\tau_{0}\mid
w_{jh},\Lambda_{-jh})$ is power law tail distributed. Hence, if any of the
scale parameters is power law tail distributed for any prior on $\beta$, then
its distribution, conditionally on $\Lambda_{-jh}$ and on the product of the
other two parameters, is power law tail distributed and, as a consequence,
$(\lambda_{jh}\mid\Lambda_{-jh})$ is power law tail distributed. Since
$f_{\lambda_{jh}\mid\Lambda_{-jh}}(\lambda)\geq c|\lambda|^{-\alpha}$ for
certain $c,\alpha$ positive constants and $|\lambda|>L$ sufficiently large, in
the same settings, we can write

$f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)=c|\lambda|^{-\alpha}\\{1+t(|\lambda|)\\},$

where $t(|\lambda|)$ is a positive function. Then,

$\displaystyle\frac{\partial[\log\\{f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)\\}]}{\partial\lambda}$
$\displaystyle=-\frac{\alpha}{\lambda}+\frac{\partial
t(\lambda)}{\partial\lambda}\qquad$ $\displaystyle\text{for}\;\lambda>L\quad$
$\displaystyle\text{and}\;L\gg 0,$
$\displaystyle\frac{\partial[\log\\{f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)\\}]}{\partial\lambda}$
$\displaystyle=\frac{\alpha}{\lambda}+\frac{\partial\\{-t(\lambda)\\}}{\partial\lambda}\qquad$
$\displaystyle\text{for}\;\lambda<-L\quad$ $\displaystyle\text{and}\;L\gg 0,$

We now consider the sign of the derivative of $t(|\lambda|)$. If
$t(|\lambda|)$ is not decreasing,

$\displaystyle\frac{\partial[\log\\{f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)\\}]}{\partial\lambda}$
$\displaystyle\geq-\frac{\alpha}{\lambda},\qquad$
$\displaystyle\text{for}\;\lambda>L\quad$ $\displaystyle\text{and}\;L\gg 0,$
$\displaystyle\frac{\partial[\log\\{f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)\\}]}{\partial\lambda}$
$\displaystyle\leq\frac{\alpha}{\lambda},\qquad$
$\displaystyle\text{for}\;\lambda<-L\quad$ $\displaystyle\text{and}\;L\gg 0,$

whereas if $t(|\lambda|)$ is decreasing, its derivative goes to zero when
$|\lambda|$ goes to infinity. Therefore,

$\displaystyle\frac{\partial[\log\\{f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)\\}]}{\partial\lambda}$
$\displaystyle\geq f_{lb}^{\prime}(\lambda)\qquad$
$\displaystyle\text{for}\;\lambda>L\quad$ $\displaystyle\text{and}\;L\gg 0,$
$\displaystyle\frac{\partial[\log\\{f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)\\}]}{\partial\lambda}$
$\displaystyle\leq-f_{lb}^{\prime}(|\lambda|)\qquad$
$\displaystyle\text{for}\;\lambda<-L\quad$ $\displaystyle\text{and}\;L\gg 0,$

where $f_{lb}^{\prime}(\lambda)<0\;\forall\,\lambda>0$ and
$\lim_{\lambda\to\infty}f_{lb}^{\prime}(|\lambda|)=0$. The proof is concluded
by using this result along with the fact that
$f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(|\lambda|)$ is decreasing when
$\lambda\rightarrow\infty$,

$\displaystyle\frac{\partial[\log\\{f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)\\}]}{\partial\lambda}$
$\displaystyle\leq 0\qquad$ $\displaystyle\text{for}\;\lambda>L\quad$
$\displaystyle\text{and}\;L\gg 0$
$\displaystyle\frac{\partial[\log\\{f_{\lambda_{jh}\mid\Lambda_{-(jh)}}(\lambda)\\}]}{\partial\lambda}$
$\displaystyle\geq 0\qquad$ $\displaystyle\text{for}\;\lambda>-L\quad$
$\displaystyle\text{and}\;L\gg 0,$

showing that the limit of the derivative for $|\lambda|\to\infty$ is equal to
zero. ∎

###### of Lemma 2.

In both the multiplicative gamma process and cumulative shrinkage process,
priors on $\Lambda$ are exchangeable within columns, that is
$\text{pr}(|\lambda_{jh}|>\epsilon\mid\gamma_{h},\tau_{0})=\zeta_{\epsilon h}$
does not depend on $j$. Then, the prior density of
$|\text{supp}_{\epsilon}(\lambda_{h})|$, conditionally on $\gamma_{h}$ and
$\tau_{0}$ is a priori distributed as a sum of independent and identically
distributed Bernoulli random variables $\text{Ber}(\zeta_{\epsilon h})$.
Furthermore, $\zeta_{\epsilon h}$ does not depend on $p$. By applying the
Chernoff’s method, we obtain

$\text{pr}\\{|\text{supp}_{\epsilon}(\lambda_{h})|<as_{p}\mid\gamma_{h},\tau_{0}\\}\leq\exp\big{\\{}ats_{p}+p\zeta_{\epsilon
h}(e^{-t}-1)\big{\\}},$

for any $t>0$ and with $1-e^{-t}>0$. Hence,

$\text{pr}\\{|\text{supp}_{\epsilon}(\lambda_{h})|>as_{p}\mid\gamma_{h},\tau_{0}\\}\geq
1-\exp[-p\\{(1-e^{-t})\zeta_{\epsilon h}-ats_{p}/p\\}],$

where the limit of the lower bound is
$\lim_{p\to\infty}1-\exp[-p\\{(1-e^{-t})\zeta_{\epsilon h}-ats_{p}/p\\}]=1$,
which concludes the proof. ∎

###### of Corollary 1.

i. It is sufficient to prove the conditions required by Theorem 1. We have
$E(\gamma_{h})=E(\vartheta_{h})E(\rho_{h})=E(\rho_{h})\,b_{\theta}/(a_{\theta}-1),$
where

$E(\rho_{h})=1-\sum_{l=1}^{h}E(w_{l})=1-\sum_{l=1}^{h-1}E(w_{l})-E(w_{h})=E(\rho_{h-1})-E(w_{h}).$

Since the random variable $w_{l}$ is obtained as a product of positive random
variables, $E(w_{l})>0$ for every $l=1,\ldots,h$. Therefore
$E(\gamma_{h})<E(\gamma_{h-1})$ for each $h=2,\ldots,H$.

ii. It is sufficient to prove the conditions required by Proposition 1. It is
straightforward to verify $E(\tau_{0})=1$ and $E(\phi_{jh})\leq 1$ for
$j=1,\ldots,p$ and $h=1,\ldots,\infty$. The column scale expectation is

$E(\gamma_{h})=E(\vartheta_{h})\left(\frac{\alpha}{1+\alpha}\right)\left(\frac{\alpha}{1+\alpha}\right)^{h-1},$

which can be written in a form $ab^{h-1}$. The elements $\sigma_{j}^{-2}$ are
gamma distributed guaranteeing finite expectation for all $j=1,\ldots,p$. ∎

###### of Corollary 2.

It is sufficient to prove the conditions required by Theorem 2. The
probability density function of the column scale $\gamma_{h}$ ($h=1,2,\ldots$)
of model (5) evaluated at a certain $\gamma>0$ is

$f_{\gamma_{h}}(\gamma)=\text{pr}(\rho_{h}=1)f_{\vartheta_{h}}(\gamma)\propto\gamma^{-a_{\theta}-1}\exp(-b_{\theta}/\gamma),$

where $f_{\vartheta_{h}}(\gamma)$ is the inverse gamma probability density
function evaluated at $\gamma$. The function
$\gamma^{-a_{\theta}-1}\exp(-b_{\theta}/\gamma)$ is of order
$O\\{\gamma^{-(a_{\theta}+1)}\\}$ as $\gamma$ goes to infinity. Since
$a_{\theta}>0$, we conclude that the column scale $\gamma_{h}$ is power law
tail distributed. The independence between $\gamma_{h}$ and $\beta_{h}$
($h=1,2,\ldots$) guarantees that the latter result hold for any possible prior
distribution $f_{\beta}$ on $\beta$. ∎

###### of Corollary 3.

It is sufficient to prove the conditions required by Theorem 3. The structured
increasing shrinkage prior is such that, for every $j=1,\ldots,p$ and $h\geq
1$, we have $g(x_{j}^{T}\beta_{h})\leq c_{p}<1$. The proof is obtained under
the assumption $c_{p}=O\\{\log(p)/p\\}$. ∎

##  Appendix S2 Simulation experiments

###  S2.1 Gibbs sampler for structured increasing shrinkage model for Gaussian
data

We can rewrite the model for $y_{ij}$ for the specific case of the structured
increasing shrinkage process and Gaussian data as

$y_{ij}=\sum_{h=1}^{\infty}\sqrt{\rho_{h}}\sqrt{\phi_{jh}}\,\lambda_{jh}^{*}\eta_{ih}+\epsilon_{ij}\qquad\lambda_{jh}^{*}\sim
N(0,\vartheta_{h}),$

where $\lambda_{jh}^{*}$ is a continuous random variable and we let
$\beta_{mh}\sim N(0,\sigma_{\beta}^{2})$. The notation $(x\mid-)$ denotes the
full conditional distribution of $x$ conditionally on everything else. Given
$H$ the number of factors of the truncated model, the sampler cycles through
the following steps.

Step S1. Update, for $i=1,\ldots,n$, the factor $\eta_{i}$ according to the
posterior full conditional

$(\eta_{i}\mid-)\sim
N_{H}\big{\\{}(I_{H}+\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}}\Sigma^{-1}\Lambda_{H})^{-1}\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}}\Sigma^{-1}y_{i},\,(I_{H}+\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}}\Sigma^{-1}\Lambda_{H})^{-1}\big{\\}}.$

Step S2. Update, for $j$ in $1,\ldots,p$, the elements of $\Sigma$, by
sampling

$(\sigma_{j}^{-2}\mid-)\sim\text{Ga}\left\\{a_{\sigma}+\frac{n}{2},\,b_{\sigma}+\frac{1}{2}\sum_{i=1}^{n}(y_{ij}-\lambda_{j}^{\mathrm{\scriptscriptstyle
T}}\eta_{i})^{2}\right\\}.$

Step S3. Update $\beta_{h}$ ($h=1,\ldots,H$) exploiting the Pólya-Gamma data-
augmentation strategy (Polson et al., 2013) and the decompostition
$\phi_{jh}=\phi_{jh}^{(L)}\phi_{jh}^{(C)}$, with
$\phi_{jh}^{(L)}\phi_{jh}^{(C)}$ independent a priori and distributed as
Ber$\\{\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle T}}\beta_{h})\\}$
and Ber$(c_{p})$, respectively.

Substep S3.1. Update $\phi_{jh}^{(L)}$, for $j=1,\ldots,p$ and $h=1,\ldots,H$,
setting $\phi_{jh}^{(L)}=1$ if $\phi_{jh}=1$ and sampling from the full
conditional distribution

$\text{pr}(\phi_{jh}^{(L)}=l)\propto\begin{cases}1-\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})\qquad\qquad\qquad\quad\,\text{for}\;l=0,\\\
\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})(1-c_{p})\qquad\,\text{for}\;l=1,\\\ \end{cases}$

if $\phi_{jh}=0$.

Substep S3.2. Let $f(y)\propto\sum_{n=0}^{\infty}(-1)^{n}A_{n}(2\pi
y^{3})^{-0.5}\exp\\{-(2n+b)^{2}(8y)^{-1}-0.5c^{2}y\\}$ indicate the
probability density function of a Pólya-Gamma distributed random variable
$y\sim\text{PG}(b,c)$. For each $h=1,\ldots,H$, generate $p$ independent
random variables $d_{j(h)}$ sampling from the full conditional distribution
$(d_{j(h)}\mid-)\sim\text{PG}(1,x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})$. Let $D_{(h)}$ denote the $p\times p$ diagonal matrix with
entries $d_{j(h)}$ ($j=1,\ldots,p$).

Substep S3.3. Define the $q\times q$ diagonal matrix
$B=\sigma^{2}_{\beta}I_{q}$. For each $h=1,\ldots,H$, update $\beta_{h}$
sampling from

$(\beta_{h}\mid-)\sim N_{q}\\{(x^{\mathrm{\scriptscriptstyle
T}}D_{(h)}x+B^{-1})^{-1}(x^{\mathrm{\scriptscriptstyle
T}}\kappa_{h}),\,(x^{\mathrm{\scriptscriptstyle T}}D_{(h)}x+B^{-1})^{-1}\\},$

where $\kappa_{h}$ is a $p$-dimensional vector with the $j$-th entry equal to
$\phi_{jh}^{(L)}-0.5$.

Step S4. Update the elements $\lambda_{jh}^{*}$ by sampling from the
independent full conditional posterior distributions of the row vectors
$\lambda_{j}^{*}=(\lambda_{j1}^{*},\ldots,\lambda_{jH}^{*})$, for
$j=1,\ldots,p$,

$(\lambda_{j}^{*}\mid-)\sim
N_{H}\big{\\{}(D^{-1}+\sigma_{j}^{-2}\eta_{(j)}^{{\mathrm{\scriptscriptstyle
T}}}\eta_{(j)})^{-1}\sigma_{j}^{-2}\eta_{(j)}^{{\mathrm{\scriptscriptstyle
T}}}y^{(j)},\,(D^{-1}+\sigma_{j}^{-2}\eta_{(j)}^{{\mathrm{\scriptscriptstyle
T}}}\eta_{(j)})^{-1}\big{\\}},$

where $\eta_{(j)}$ is the $n\times H$ matrix such that the generic element is
$\eta_{(j)ih}=\eta_{ih}\sqrt{\rho_{h}}\sqrt{\phi_{jh}}$,
$D^{-1}=\text{diag}(\vartheta_{1}^{-1},\ldots,\vartheta_{H}^{-1})$ and
$y^{(j)}=(y_{1j},\ldots,y_{nj})^{\mathrm{\scriptscriptstyle T}}$. Set
$\lambda_{jh}=\lambda_{jh}^{*}\sqrt{\rho_{h}}\sqrt{\phi_{jh}}$.

Step S5. Update the column scales $\gamma_{h}$ (for $h=1,\ldots,H$), following
the substeps below and setting $\gamma_{h}=\vartheta_{h}\rho_{h}$.
Consistently with Legramanti et al. (2020), define the independent indicators
$z_{h}$ ($h=1,\ldots,p$) with prior $\text{pr}(z_{h}=l)=w_{l}$.

Substep S5.1. Update the augmented data $z_{h}$ by sampling from the full
conditional distribution

$\text{pr}(z_{h}=l)\propto\begin{cases}w_{l}\,\prod_{i=1}^{n}\prod_{j=1}^{p}N(y_{ij};\mu_{ijh}^{(0)},\sigma_{j}^{2})\qquad\qquad\text{for}\quad
l=1,\ldots,h\\\
w_{l}\,\prod_{i=1}^{n}\prod_{j=1}^{p}N(y_{ij};\mu_{ijh}^{(1)},\sigma_{j}^{2})\qquad\qquad\text{for}\quad
l=h+1,\ldots,H,\\\ \end{cases}$ (10)

where $N(x;\mu,\sigma^{2})$ indicates the Gaussian probability density
function with mean $\mu$ and variance $\sigma^{2}$. The mean values
$\mu_{ijh}^{(0)}$ and $\mu_{ijh}^{(1)}$ are defined according to
$\mu_{ijh}^{(z)}=\sum_{l\neq
h}^{H}\sqrt{\rho_{l}}\,\sqrt{\phi_{jl}}\lambda_{jl}^{*}\eta_{il}+\sqrt{z}\,\sqrt{\phi_{jh}}\lambda_{jh}^{*}\eta_{ih}$.
Set $\rho_{h}=1$ if $z_{h}>h$, else $\rho_{h}=0$.

Substep S5.2. For $h=1,\ldots,H$, update $\vartheta_{h}^{-1}$ sampling from
$\text{Ga}(a_{\theta}+0.5p,b_{\theta}+0.5\sum_{j=1}^{p}\lambda_{jh}^{*\,2})$.

Substep S5.3. For $l=1,\ldots,H-1$, sample $v_{l}$ from

$(v_{l}\mid-)\sim\text{Be}\big{\\{}1+\sum_{h=1}^{H}\mathbbm{1}{(z_{h}=l)},\alpha+\mathbbm{1}{(z_{h}>l)}\big{\\}},$

set $v_{H}=1$ and update $w_{l}=v_{l}\prod_{m=1}^{l-1}(1-v_{m})$, for
$l=1,\ldots,H$.

Step S6. Update independently the local scales, for $j=1,\ldots,p$ and
$h=1,\ldots,H$, by sampling from the full conditional distributions

$\text{pr}(\phi_{jh}=u)\propto\begin{cases}\\{1-\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})\,c_{p}\\}\,\prod_{i=1}^{n}N(y_{ij};\mu_{ijh}^{(u)},\sigma_{j}^{2})\quad\;\text{for}\;u=0\\\
\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})\,c_{p}\prod_{i=1}^{n}N(y_{ij};\mu_{ijh}^{(u)},\sigma_{j}^{2})\qquad\qquad\text{for}\;u=1.\end{cases}$

with $\mu_{ijh}^{(u)}=\sum_{l\neq
h}^{H}\sqrt{\rho_{l}}\,\sqrt{\phi_{jl}}\lambda_{jl}^{*}\eta_{il}+\sqrt{\rho_{h}}\,\sqrt{u}\lambda_{jh}^{*}\eta_{ih}$.

###  S2.2 Simulation settings

The results reported in Section 4 are obtained running the algorithms for
25000 iterations discarding the first 10000 iterations. Then, we thin the
Markov chain, saving every $5$-th sample. We adapt the number of active
factors at iteration $t$ with probability $p(t)=\exp(-1-5\,10^{-4}t)$. We set
$a_{\sigma}=1$ and $b_{\sigma}=0.3$. In the structured increasing shrinkage
algorithm, we choose the offset constant $c_{p}=2e\log(p)/p$ which belongs to
$(0,1)$ for every $p\geq 15$.

In scenario $d$, the meta covariates in matrix $x_{0}$ are a categorical
variable with four balanced categories, a continuous variable sampled from a
multivariate Gaussian distribution, and a continuous variable where the $p$
elements are sampled from $p$ gamma distributions.

To infer the structural zeros within each column of $\Lambda$ in the
cumulative shrinkage process and in the multiplicative gamma process, we set
$\lambda_{jh}$ to zero when $|\lambda_{jh}|$ ($j=1,\ldots,p$) is under a
certain threshold. We choose the threshold equal to 0.05, which is consistent
with the value of the hyperparameter $\theta_{\infty}$ used in the cumulative
shrinkage process.

To address column order ambiguity and label switching, we compute the mean
classification error only after having ordered the columns of $\Lambda^{(t)}$
(for $t=1,\ldots,S$), for each model, increasingly with respect to the number
of zero entries identified.

###  S2.3 Simulation results

We report additional results for the simulation study of Section 4 of the main
paper.

Table S3: Median and interquartile range of the LPML, Cov. MSE and of
$E(H_{a}\mid y)$ computed in 25 replications assuming Scenario b and several
combinations of $(p,k,s)$

| $(p,k,s)$ | MGP | CUSP | SIS
---|---|---|---|---
|  | Q0.5 | IQR | Q0.5 | IQR | Q0.5 | IQR
LPML |  00(16,4,0.6) |  0-28.20 | 0.33 |  0-28.20 | 0.33 |  0-28.17 | 0.32
|  00(32,8,0.4) |  0-56.95 | 0.53 |  0-57.00 | 0.51 |  0-56.80 | 0.49
|  0(64,12,0.3) | -111.35 | 0.70 | -111.71 | 0.74 | -110.76 | 0.89
| (128,16,0.2) | -211.65 | 0.74 | -215.94 | 1.57 | -210.19 | 0.86
Cov. MSE |  00(16,4,0.6) |  000.25 | 0.12 |  000.25 | 0.12 |  000.23 | 0.10
|  00(32,8,0.4) |  000.32 | 0.08 |  000.33 | 0.10 |  000.30 | 0.12
|  0(64,12,0.3) |  000.37 | 0.10 |  000.43 | 0.11 |  000.22 | 0.09
| (128,16,0.2) |  000.23 | 0.03 |  000.32 | 0.04 |  000.09 | 0.01
$E(H_{a}\mid y)$ |  00(16,4,0.6) |  008.91 | 1.52 |  004.00 | 0.00 |  004.00 | 0.00
|  00(32,8,0.4) |  011.27 | 1.48 |  007.00 | 1.00 |  008.00 | 0.00
|  0(64,12,0.3) |  014.72 | 1.49 |  011 .00 | 0.00 |  012.00 | 0.00
| (128,16,0.2) |  017.16 | 0.81 |  012.00 | 1.75 |  016.00 | 0.00

LPML, logarithm of the pseudo-marginal likelihood; Cov. MSE, covariance mean
squared error; CUSP, cumulative increasing shrinkage process; MGP,
multiplicative gamma process; SIS, structured increasing shrinkage process;
Q0.5, median; IQR, interquartile range.

In scenario b, we also apply the method proposed by Ročková and George (2016),
which is referred to as parameter expanded likelihood expected maximization.
Hyperparameters are set as suggested by the authors. This approach focuses on
finding a sparse mode based on an over-parameterized factor model. The
performance in terms of mean squared error in covariance estimation and
classification error in detecting sparsity in $\Lambda$ is reported in Table
S2. The results are not competitive with the other approaches we have
considered.

Table S4: Median and interquartile range of the Cov. MSE and CE computed
applying the parameter expanded likelihood expected maximization method of
Ročková and George (2016) in 25 replications under Scenario b and several
combinations of $(p,k,s)$

$(p,k,s)$ | Cov. MSE | CE
---|---|---
| Q0.5 | IQR | Q0.5 | IQR
00(16,4,0.6) | 0.55 | 0.16 | 0.77 | 0.14
00(32,8,0.4) | 0.66 | 0.11 | 0.76 | 0.10
0(64,12,0.3) | 0.64 | 0.11 | 0.90 | 0.10
(128,16,0.2) | 0.42 | 0.03 | 1.06 | 0.20

Cov. MSE, covariance mean squared error; CE, classification error; PXLEM,
parameter expanded likelihood expected maximization; Q0.5, median; IQR,
interquartile range.

Table S5: Median and interquartile range of the LPML, Cov. MSE, $E(H_{a}\mid
y)$ and MCE computed in 25 replications assuming Scenario c and several
combinations of $(p,k,s)$

| $(p,k,s)$ | MGP | CUSP | SIS
---|---|---|---|---
|  | Q0.5 | IQR | Q0.5 | IQR | Q0.5 | IQR
LPML |  00(16,4,0.6) |  0-27.62 | 0.24 |  0-27.62 | 0.25 |  0-27.59 | 0.24
|  00(32,8,0.4) |  0-56.16 | 0.64 |  0-56.22 | 0.51 |  0-55.89 | 0.59
|  0(64,12,0.3) | -109.64 | 0.69 | -110.67 | 0.71 | -109.06 | 0.65
| (128,16,0.2) | -209.57 | 0.88 | -214.19 | 1.76 | -208.34 | 1.04
Cov. MSE |  00(16,4,0.6) |  000.30 | 0.10 |  000.29 | 0.09 |  000.26 | 0.11
|  00(32,8,0.4) |  000.77 | 0.26 |  000.72 | 0.18 |  000.80 | 0.43
|  0(64,12,0.3) |  001.01 | 0.35 |  000.94 | 0.21 |  001.20 | 1.22
| (128,16,0.2) |  000.78 | 0.18 |  000.87 | 0.21 |  000.35 | 0.48
$E(H_{a}\mid y)$ |  00(16,4,0.6) |  008.38 | 1.80 |  003.44 | 1.00 |  004.00 | 0.00
|  00(32,8,0.4) |  010.38 | 1.12 |  005.05 | 0.91 |  008.00 | 1.00
|  0(64,12,0.3) |  013.67 | 1.20 |  008.00 | 0.92 |  012.00 | 0.00
| (128,16,0.2) |  016.56 | 0.83 |  009.00 | 0.00 |  016.00 | 0.00
MCE |  00(16,4,0.6) |  000.98 | 0.17 |  000.53 | 0.20 |  000.24 | 0.06
|  00(32,8,0.4) |  000.65 | 0.07 |  000.44 | 0.08 |  000.19 | 0.07
|  0(64,12,0.3) |  000.59 | 0.04 |  000.48 | 0.04 |  000.18 | 0.06
| (128,16,0.2) |  000.48 | 0.02 |  000.44 | 0.01 |  000.06 | 0.10

LPML, logarithm of the pseudo-marginal likelihood; Cov. MSE, covariance mean
squared error; MCE, mean classification error; CUSP, cumulative increasing
shrinkage process; MGP, multiplicative gamma process; SIS, structured
increasing shrinkage process; Q0.5, median; IQR, interquartile range.

Table S6: Median and interquartile range of the LPML, Cov. MSE, $E(H_{a}\mid
y)$ and MCE computed in 25 replications assuming Scenario d and several
combinations of $(p,k,s)$

| $(p,k,s)$ | MGP | CUSP | SIS | SISmc
---|---|---|---|---|---
|  | Q0.5 | IQR | Q0.5 | IQR | Q0⋅5 | IQR | Q0.5 | IQR
LPML |  00(16,4,0.6) |  0-27.74 | 0.43 |  0-27.75 | 0.43 |  0-27.71 | 0.40 |  0-27.73 | 0.40
|  00(32,8,0.4) |  0-56.25 | 0.69 |  0-56.35 | 0.72 |  0-56.16 | 0.68 |  0-56.12 | 0.65
|  0(64,12,0.3) | -109.72 | 0.61 | -110.54 | 0.88 | -109.27 | 0.46 | -109.16 | 0.57
| (128,16,0.2) | -209.60 | 0.48 | -213.50 | 1.21 | -208.11 | 0.42 | -208.03 | 0.47
Cov. MSE |  00(16,4,0.6) |  000.31 | 0.11 |  000.30 | 0.14 |  000.28 | 0.14 |  000.27 | 0.16
|  00(32,8,0.4) |  000.70 | 0.25 |  000.71 | 0.18 |  000.75 | 0.22 |  000.79 | 0.78
|  0(64,12,0.3) |  001.03 | 0.43 |  000.91 | 0.29 |  001.51 | 0.59 |  001.16 | 0.84
| (128,16,0.2) |  000.93 | 0.49 |  000.90 | 0.33 |  001.49 | 1.21 |  001.28 | 1.81
$E(H_{a}\mid y)$ |  00(16,4,0.6) |  008.60 | 0.64 |  003.96 | 0.80 |  004.00 | 0.00 |  004.00 | 0.00
|  00(32,8,0.4) |  010.71 | 1.24 |  005.75 | 1.00 |  007.00 | 1.00 |  008.00 | 0.00
|  0(64,12,0.3) |  013.93 | 1.37 |  008.00 | 0.92 |  012.00 | 0.00 |  012.00 | 0.00
| (128,16,0.2) |  016.56 | 0.88 |  009.00 | 0.00 |  016.00 | 1.00 |  016.00 | 1.00
MCE |  00(16,4,0.6) |  000.94 | 0.13 |  000.64 | 0.19 |  000.26 | 0.08 |  000.23 | 0.13
|  00(32,8,0.4) |  000.67 | 0.10 |  000.49 | 0.09 |  000.20 | 0.08 |  000.20 | 0.10
|  0(64,12,0.3) |  000.58 | 0.05 |  000.47 | 0.04 |  000.21 | 0.06 |  000.21 | 0.08
| (128,16,0.2) |  000.49 | 0.02 |  000.43 | 0.02 |  000.18 | 0.11 |  000.17 | 0.11

LPML, logarithm of the pseudo-marginal likelihood; Cov. MSE, covariance mean
squared error; MCE, mean classification error; CUSP, cumulative increasing
shrinkage process; MGP, multiplicative gamma process; SIS, structured
increasing shrinkage process; SISmc, structured increasing shrinkage process
with meta covariates; Q0.5, median; IQR, interquartile range.

Figure S6: Boxplots of the logarithm of the pseudo-marginal likelihood for all
combinations $(p,k)$ in scenario a (top left panel) and scenario b (top right
panel), of the covariance mean square error in scenario a (bottom left panel),
and of the mean classification error in scenario b (bottom right panel). LPML,
logarithm of the pseudo-marginal likelihood; Cov. MSE, covariance mean squared
error; MCE, mean classification error; CUSP, cumulative increasing shrinkage
process; MGP, multiplicative gamma process; SIS, structured increasing
shrinkage process.

Figure S7: Boxplots of the logarithm of the pseudo-marginal likelihood and of
the mean classification error of each model for all combinations of $(p,k,s)$
in Scenario c (left panel) and Scenario d (right panel). LPML, logarithm of
the pseudo-marginal likelihood; MCE, mean classification error; CUSP,
cumulative shrinkage process; MGP, multiplicative gamma process; SIS,
structured increasing shrinkage process; SIS mc, structured increasing
shrinkage process with meta covariates.

###  S2.4 Simulation study of sensitivity to hyperparameters and truncation
level

We conduct further simulation experiments to assess the impact of some
hyperparameters on relevant prior and posterior summaries. Figure S8 displays
the prior distribution, obtained simulating 10,000 samples from the prior, of
the proportion of variance explained by the structured increasing shrinkage
factor model for varying $\alpha$,
$\\{E(\sigma^{-2}),\text{var}(\sigma^{-2})\\}$, and
$\\{E(\vartheta_{h}^{-1}),\text{var}(\vartheta_{h}^{-1})\\}$. The
hyperparameter $\alpha$, representing the expected number of factors,
positively affects the proportion of explained variance. The influence of the
hyperparameters regulating the distribution of $\vartheta_{h}$ is even
clearer, with concentrated prior on a large value of $E(\vartheta^{-1})$,
inducing a smaller proportion of variance explained by the factor model. The
role of $\\{E(\sigma^{-2}),\text{var}(\sigma^{-2})\\}$ is less clear, but
suggests that sufficiently large mean and variance can guarantee higher
flexibility.

Figure S8: Boxplots of the prior proportion of variance explained by the
factor model $\text{tr}(\Lambda\Lambda^{\mathrm{\scriptscriptstyle
T}})/\text{tr}(\Omega)$. The quantity is obtained simulating 10,000 samples
from the prior distribution with varying values of the parameters. The
horizontal axis characterize the effect of
$\\{E(\sigma^{-2}),\text{var}(\sigma^{-2})\\}$; differences for
$\alpha\in(5,10,20)$ are reported in each column; differences for
$\\{E(\vartheta^{-1}),\text{var}(\vartheta^{-1})\\}\in\\{(2,2),(1,0.5),(0.5,0.125)\\}$
are reported in each row.

The latter comment is confirmed by the study of the impact of $\alpha$ and
$\\{E(\sigma^{-2}),\text{var}(\sigma^{-2})\\}$ on the posterior bound of the
truncation error and on the posterior distribution of the proportion of
variance explained by the factor model. Specifically, we generate synthetic
data sets with $n=100$ observations with dimension $p=50$ from the Gaussian
linear factor model $y_{i}\sim N_{p}(0,\Lambda_{0}\Lambda_{0}^{T}+I_{p})$,
with $\Lambda_{0}$ a sparse $p\times k$ matrix with $k=50$. We randomly set
two thirds of the elements of $\Lambda_{0}$ equal to zero, drawing the non
zero elements from a Gaussian distribution with mean zero and variances
$\theta_{h}$ sampled from an inverse gamma distribution
$\theta_{h}^{-1}\sim\text{Ga}(2,2)$. We keep the number of active factors $H$
fixed at $50$, and set $a_{\theta}=b_{\theta}=2$ and $c_{p}=2e\log(p)/p$. We
run the Gibbs algorithm for the structured increasing shrinkage model for
15000 iterations, discarding the first 5000 iterations. Then, we thin the
Markov chain, saving every 5-th sample.

In Figure S9 the sampled posterior distribution of the proportion of variance
explained by the factor model
$\text{tr}(\Lambda\Lambda^{\mathrm{\scriptscriptstyle T}})/\text{tr}(\Omega)$
is reported for varying $\alpha$ and
$\\{E(\sigma^{-2}),\text{var}(\sigma^{-2})\\}$. The same proportion computed
on the matrices generating the data is
$\text{tr}(\Lambda_{0}\Lambda_{0}^{\mathrm{\scriptscriptstyle
T}})/\text{tr}(\Omega_{0})=0.966$, with
$\Omega_{0}=\Lambda_{0}\Lambda_{0}^{\mathrm{\scriptscriptstyle T}}+I_{50}$. A
sufficiently concentrated prior on a large value of $E(\sigma^{-2})$ seems
more suitable to model such data, even if we have incorrect expectations on
the number of factors, i.e. $\alpha$ set small.

Figure S9: Boxplots representing the simulated posterior distribution of the
proportion of variance explained by the factor model
$\text{tr}(\Lambda\Lambda^{\mathrm{\scriptscriptstyle T}})/\text{tr}(\Omega)$
for varying $\alpha$ and $\\{E(\sigma^{-2}),\text{var}(\sigma^{-2})\\}$. The
dashed lines represent the proportion computed on the true value of $\Lambda$
and $\Omega$.

Figure S10 displays the Monte Carlo approximation of the posterior probability
of truncation error $\text{pr}\\{\text{tr}(\Omega_{H})/\text{tr}(\Omega)<T\\}$
for different values of $H$ and $T$ and varying $\alpha$ and
$\\{E(\sigma^{-2}),\text{var}(\sigma^{-2})\\}$. If $\Lambda_{0}$ is sparse, a
small value of $\alpha$ induces good approximations even with $H$ smaller than
the true number of factors. The inferred sparsity pattern in $\Lambda$ is
robust to the prior distribution for $\sigma^{-2}$.

Figure S10: Monte Carlo approximation of the posterior probability of
truncation error $\text{pr}(r<T\mid y)$, with
$r=\text{tr}(\Omega_{H})/\text{tr}(\Omega)$, at varying of $H$. The quantity
is computed for $T$ equal to $0.75$ (—), $0.9$ (- - -), and $0.95$ (– – –) and
varying $\alpha\in(5,10,20)$ over the columns and
$\\{E(\sigma^{-2}),\text{var}(\sigma^{-2})\\}\in\\{(3.33,11.11),(1,1),(0.33,0.11)\\}$
over the rows of the figure.

##  Appendix S3 Finnish bird co-occurrence application

###  S3.1 Gibbs algorithms of probit structured increasing shrinkage model

In case of probit data (see Section 5 of the main paper) and the structured
increasing shrinkage process, we can rewrite the latent model for $z_{ij}$ as

$\displaystyle z_{ij}$ $\displaystyle=w_{i}^{T}\mu_{j}+\epsilon_{ij},$
$\displaystyle\epsilon_{ij}$
$\displaystyle=\sum_{h=1}^{\infty}\sqrt{\rho_{h}}\sqrt{\phi_{jh}}\,\lambda_{jh}^{*}\eta_{ih}+\varepsilon_{ij},\qquad\lambda_{jh}^{*}\sim
N(0,\vartheta_{h}),\qquad\varepsilon_{ij}\sim N(0,1),$

where $\lambda_{jh}^{*}$ is an absolutely continuous random variable. Let the
notation $(x\mid-)$ denote the full conditional distribution of $x$
conditionally on everything else. Given $H$ the number of factors of the
truncated model, the sampler cycles through the following steps.

Step S1. Update $\mu_{j}$, for every $j=1,\ldots,p,$ by sampling from the
independent full conditional posterior distributions

$(\mu_{j}\mid-)\sim
N_{c}\big{[}(\sigma_{\mu}^{-2}I_{c}+w^{\mathrm{\scriptscriptstyle
T}}w)^{-1}\\{w^{\mathrm{\scriptscriptstyle
T}}(z^{(j)}-\eta\lambda_{j})+bx_{j}\\},(\sigma_{\mu}^{-2}I_{c}+w^{\mathrm{\scriptscriptstyle
T}}w)^{-1}\big{]},$

where $z^{(j)}=(z_{1j},\ldots,z_{nj})^{\mathrm{\scriptscriptstyle T}}$ and
$\eta=(\eta_{1},\ldots,\eta_{n})^{\mathrm{\scriptscriptstyle T}}$.

Step S2. Update $b_{l}$ ($l=1,\ldots,c$) sampling from conditionally
independent posteriors

$(b_{l}\mid-)\sim
N_{q}\big{\\{}(\sigma_{b}^{-2}I_{q}+\sigma_{\mu}^{-2}x^{\mathrm{\scriptscriptstyle
T}}x)^{-1}\sigma_{\mu}^{-2}(x^{\mathrm{\scriptscriptstyle
T}}\mu^{(l)}),\,(\sigma_{b}^{-2}I_{q}+\sigma_{\mu}^{-2}x^{\mathrm{\scriptscriptstyle
T}}x)^{-1}\big{\\}},$

where $\mu^{(l)}=(\mu_{1l},\ldots,\mu_{pl})^{\mathrm{\scriptscriptstyle T}}$

Step S3. Update the elements $z_{ij}$ ($i=1,\ldots,n$; $j=1\ldots,p$) sampling
independently from the truncated normal

$(z_{ij}\mid-)\sim TN(\lambda_{j}^{\mathrm{\scriptscriptstyle
T}}\eta_{i}+w_{i}^{T}\mu_{j},1,l_{ij},u_{ij}),$

where the lower bound $l_{ij}$ is equal to $0$ if $y_{ij}=1$ and $-\infty$
otherwise. The upper bound $u_{ij}=0$ if $y_{ij}=0$ and $\infty$ otherwise.
Then, set $\epsilon=z-w\mu$.

Step S4. Update, for $i=1,\ldots,n$, the factor $\eta_{i}$ according to the
posterior full conditional

$(\eta_{i}\mid-)\sim
N_{H}\big{\\{}(I_{H}+\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}}\Lambda_{H})^{-1}\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}}\epsilon_{i},\,(I_{H}+\Lambda_{H}^{\mathrm{\scriptscriptstyle
T}}\Lambda_{H})^{-1}\big{\\}}.$

Step S5. Update $\beta_{h}$ ($h=1,\ldots,H$) exploiting the Pólya-Gamma data-
augmentation strategy (Polson et al., 2013) and the decompostition
$\phi_{jh}=\phi_{jh}^{(L)}\phi_{jh}^{(C)}$, with
$\phi_{jh}^{(L)}\phi_{jh}^{(C)}$ independent a priori and distributed as
Ber$\\{\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle T}}\beta_{h})\\}$
and Ber$(c_{p})$, respectively.

Substep S5.1. Update $\phi_{jh}^{(L)}$, for $j=1,\ldots,p$ and $h=1,\ldots,H$,
setting $\phi_{jh}^{(L)}=1$ if $\phi_{jh}=1$ and sampling from the full
conditional distribution

$\text{pr}(\phi_{jh}^{(L)}=l)\propto\begin{cases}1-\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})\qquad\qquad\qquad\quad\,\text{for}\;l=0,\\\
\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})(1-c_{p})\qquad\,\text{for}\;l=1,\\\ \end{cases}$

if $\phi_{jh}=0$.

Substep S5.2. Let $f(y)\propto\sum_{n=0}^{\infty}(-1)^{n}A_{n}(2\pi
y^{3})^{-0.5}\exp\\{-(2n+b)^{2}(8y)^{-1}-0.5c^{2}y\\}$ indicate the
probability density function of a Pólya-Gamma distributed random variable
$y\sim\text{PG}(b,c)$. For each $h=1,\ldots,H$, generate $p$ independent
random variables $d_{j(h)}$ sampling from the full conditional distribution
$(d_{j(h)}\mid-)\sim\text{PG}(1,x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})$. Let $D_{(h)}$ denote the $p\times p$ diagonal matrix with
entries $d_{j(h)}$ ($j=1,\ldots,p$).

Substep S5.3. Define the $q\times q$ diagonal matrix
$B=\sigma^{2}_{\beta}I_{q}$. For each $h=1,\ldots,H$, update $\beta_{h}$
sampling from

$(\beta_{h}\mid-)\sim N_{q}\\{(x^{\mathrm{\scriptscriptstyle
T}}D_{(h)}x+B^{-1})^{-1}(x^{\mathrm{\scriptscriptstyle
T}}\kappa_{h}),\,(x^{\mathrm{\scriptscriptstyle T}}D_{(h)}x+B^{-1})^{-1}\\},$

where $\kappa_{h}$ is the $p$-dimensional vector with the $j$-th entry equal
to $\phi_{jh}^{(L)}-0.5$.

Step S6. Update the elements $\lambda_{jh}^{*}$ by sampling from the
independent full conditional posterior distributions of the rows vector
$\lambda_{j}^{*}=(\lambda_{j1}^{*},\ldots,\lambda_{jH}^{*})$, for
$j=1,\ldots,p$,

$(\lambda_{j}^{*}\mid-)\sim
N_{H}\big{\\{}(D^{-1}+\eta_{(j)}^{{\mathrm{\scriptscriptstyle
T}}}\eta_{(j)})^{-1}\eta_{(j)}^{{\mathrm{\scriptscriptstyle
T}}}\epsilon^{(j)},\,(D^{-1}+\eta_{(j)}^{{\mathrm{\scriptscriptstyle
T}}}\eta_{(j)})^{-1}\big{\\}},$

where $\eta_{(j)}$ is the $n\times H$ matrix such that the generic element is
$\eta_{(j)ih}=\eta_{ih}\sqrt{\rho_{h}}\sqrt{\phi_{jh}}$,
$D^{-1}=\text{diag}(\vartheta_{1}^{-1},\ldots,\vartheta_{H}^{-1})$ and
$\epsilon^{(j)}=(\epsilon_{1j},\ldots,\epsilon_{nj})^{\mathrm{\scriptscriptstyle
T}}$. Set $\lambda_{jh}=\lambda_{jh}^{*}\sqrt{\rho_{h}}\sqrt{\phi_{jh}}$.

Step S7. Update the column scales $\gamma_{h}$ (for $h=1,\ldots,H$), following
the substeps below and setting $\gamma_{h}=\vartheta_{h}\rho_{h}$.
Consistently with Legramanti et al. (2020), define the independent indicators
$u_{h}$ ($h=1,\ldots,p$) with prior $\text{pr}(u_{h}=l)=w_{l}$.

Substep S7.1. Update the augmented data $u_{h}$ by sequentially sampling from
the full conditional distribution

$\text{pr}(u_{h}=l)\propto\begin{cases}w_{l}\,\prod_{i=1}^{n}\prod_{j=1}^{p}N(\epsilon_{ij};\mu_{ijh}^{(0)},\sigma_{j}^{2})\qquad\qquad\text{for}\quad
l=1,\ldots,h\\\
w_{l}\,\prod_{i=1}^{n}\prod_{j=1}^{p}N(\epsilon_{ij};\mu_{ijh}^{(1)},\sigma_{j}^{2})\qquad\qquad\text{for}\quad
l=h+1,\ldots,H.\\\ \end{cases}$ (11)

The mean values $\mu_{ijh}^{(0)}$ and $\mu_{ijh}^{(1)}$ are defined according
to $\mu_{ijh}^{(u)}=\sum_{l\neq
h}^{H}\sqrt{\rho_{l}}\,\sqrt{\phi_{jl}}\lambda_{jl}^{*}\eta_{il}+\sqrt{u}\,\sqrt{\phi_{jh}}\lambda_{jh}^{*}\eta_{ih}$.
Set $\rho_{h}=1$ if $u_{h}>h$, else $\rho_{h}=0$.

Substep S7.2. For $h=1,\ldots,H$, update $\vartheta_{h}^{-1}$ sampling from
$\text{Ga}(a_{\theta}+0.5p,b_{\theta}+0.5\sum_{j=1}^{p}\lambda_{jh}^{*\,2})$.

Substep S7.3. For $l=1,\ldots,H-1$, sample $v_{l}$ from

$(v_{l}\mid-)\sim\text{Be}\big{\\{}1+\sum_{h=1}^{H}\mathbbm{1}{(u_{h}=l)},\alpha+\mathbbm{1}{(u_{h}>l)}\big{\\}},$

set $v_{H}=1$ and update $w_{l}=v_{l}\prod_{m=1}^{l-1}(1-v_{m})$, for
$l=1,\ldots,H$.

Step S8. Update the local scales, independently for $j=1,\ldots,p$ and
sequentially over $h=1,\ldots,H$, by sampling from the full conditional
distributions

$\text{pr}(\phi_{jh}=u)\propto\begin{cases}\\{1-\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})\,c_{p}\\}\,\prod_{i=1}^{n}N(\epsilon_{ij};\mu_{ijh}^{(u)},1)\quad\;\text{for}\;u=0\\\
\text{logit}^{-1}(x_{j}^{\mathrm{\scriptscriptstyle
T}}\beta_{h})\,c_{p}\prod_{i=1}^{n}N(\epsilon_{ij};\mu_{ijh}^{(u)},1)\qquad\qquad\text{for}\;u=1.\end{cases}$

with $\mu_{ijh}^{(u)}=\sum_{l\neq
h}^{H}\sqrt{\rho_{l}}\,\sqrt{\phi_{jl}}\lambda_{jl}^{*}\eta_{il}+\sqrt{\rho_{h}}\,\sqrt{u}\lambda_{jh}^{*}\eta_{ih}$.

The results reported in Section 5 are obtained running the algorithm for 40000
iterations discarding the first 20000 iterations. Then, we thin the Markov
Chain, saving every $5$-th sample. We adapt the number of active factors at
iteration $t$ with probability $p(t)=\exp(-1-2.5\,10^{(-4)}t)$ and, given the
high value of $p$ considered, we choose the offset constant
$c_{p}=2e\log(p)/p$ which belongs to $(0,1)$ for every $p\geq 15$.

###  S3.2 Gibbs chains mixing

Figure S11: Chain plots of the marginal posterior samples of 12 mean
coefficients of the matrix $\mu$ obtained by the Gibbs sampler, discarding the
first 20000 iterations and saving every $5$-th sample. Figure S12: Chain plots
of the marginal posterior samples of six elements of the covariance matrix
obtained by the Gibbs sampler, discarding the first 20000 iterations and
saving every $5$-th sample.

###  S3.3 Figures

Figure S13: Posterior mean of $\mu$ and $b$ for the structured increasing
shrinkage model; rows of left matrix refer to the 50 birds species, and rows
of right matrix to the ten species traits. Broadle: broadleaved forests;
Conifer: coniferous forests; Temp: temperature.

Figure S14: Maps of the sampling units in Finland coloured accordingly to the
values of the first and the third latent factors sampled at iteration $t^{*}$.
Red and blue spots represent the environments with positive and negative
values of the factors, respectively.