Order-invariant prior specification in Bayesian factor analysis
In (exploratory) factor analysis, the loading matrix is identified only up to orthogonal rotation. For identifiability, one thus often takes the loading matrix to be lower triangular with positive diagonal entries. In Bayesian inference, a standard practice is then to specify a prior under which the loadings are independent, the off-diagonal loadings are normally distributed, and the diagonal loadings follow a truncated normal distribution. This prior specification, however, depends in an important way on how the variables and associated rows of the loading matrix are ordered. We show how a minor modification of the approach allows one to compute with the identifiable lower triangular loading matrix but maintain invariance properties under reordering of the variables.
factanal_v5.tex
\documentclass[10pt]{amsart}
 
\usepackage{amsthm,amsmath,amssymb}
\usepackage{graphicx}
\usepackage[comma, sort&compress]{natbib}
\def\bibfont{\small}
\usepackage{paralist}
%% \graphicspath{}
\RequirePackage{hyperref}


%% \usepackage[left=1.5in,right=1.5in,top=1.5in,bottom=1.5in]{geometry} 

\newtheorem{theorem}{Theorem}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{lemma}[theorem]{Lemma} 
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{question}[theorem]{Question}
\newtheorem{problem}[theorem]{Problem}
\theoremstyle{definition} 
\newtheorem{definition}{Definition}
\theoremstyle{remark} 
\newtheorem{example}{Example}
\newtheorem{remark}{Remark}

%% \newtheorem{theorem}{Theorem}
%% \newtheorem{lemma}{Lemma}
%% \newtheorem{proposition}{Proposition}
%% \newtheorem{corollary}{Corollary}
%% \newtheorem{question}{Question}
%% \newtheorem{problem}{Problem}
%% \newtheorem{conjecture}{Conjecture}
%% \theoremstyle{definition}
%% \newtheorem{definition}{Definition}
%% \newtheorem{example}{Example}
%% \newtheorem*{remark}{Remark}
%% \newtheorem*{trekrule}{Trek Rule}


\numberwithin{equation}{section}
\numberwithin{theorem}{section}
%% \numberwithin{lemma}{section}
%% \numberwithin{proposition}{section}
%% \numberwithin{corollary}{section}
\numberwithin{example}{section}
\numberwithin{definition}{section}
%% \numberwithin{remark}{section}

\numberwithin{figure}{section}
\renewcommand{\theenumi}{\roman{enumi}}

\newcommand{\an}{{\rm An}}       %for ancestors
\newcommand{\de}{{\rm de}}       %for descendants
\newcommand{\dis}{{\rm dis}}     % for district
\newcommand{\ch}{{\rm ch}}       %  for children
\newcommand{\pa}{{\rm pa}}       % for parents
\newcommand{\nei}{{\rm ne}}      % for neighbours NOTE   \nei not \ne
\newcommand{\mb}{{\rm mb}}       %for Markov blanket
\newcommand{\sib}{{\rm sib}}       %for siblings
\newcommand{\htr}{{\rm htr}}       %for half-trek-reachable nodes (dependencies)

\newcommand{\lhs}[1]{\mathrm{Left}\left({#1}\right)}
\newcommand{\rhs}[1]{\mathrm{Right}\left({#1}\right)}
\newcommand{\tp}[1]{\mathrm{Top}\left({#1}\right)}
\newcommand{\tops}{\mathrm{Tops}}
\newcommand{\tpd}[1]{\mathrm{Top_d}\left({#1}\right)}
\newcommand{\bi}{\leftrightarrow}

\DeclareMathOperator{\diag}{diag}



%----------------------------------------------------------
% Reference macros
%----------------------------------------------------------
\newcommand{\rmkref}[1]{Remark~\ref{rmk:#1}}
\newcommand{\assumpref}[1]{Assumption~\ref{assump:#1}}
\newcommand{\assumpsref}[1]{Assumptions~\ref{assump:#1}}
\newcommand{\assumpssref}[1]{\ref{assump:#1}}
\newcommand{\figref}[1]{Figure~\ref{fig:#1}}
\newcommand{\figsref}[1]{Figures~\ref{fig:#1}}
\newcommand{\figssref}[1]{\ref{fig:#1}}
\newcommand{\chapref}[1]{Chapter~\ref{chap:#1}}
\newcommand{\secref}[1]{Section~\ref{sec:#1}}
\newcommand{\appref}[1]{Appendix~\ref{app:#1}}
\newcommand{\eqssref}[1]{\ref{eq:#1}}
\newcommand{\defref}[1]{Definition~\ref{def:#1}}
\newcommand{\exref}[1]{Example~\ref{ex:#1}}
\newcommand{\clmref}[1]{Claim~\ref{clm:#1}}
\newcommand{\lemref}[1]{Lemma~\ref{lem:#1}}
\newcommand{\lemsref}[1]{Lemmas~\ref{lem:#1}}
\newcommand{\lemssref}[1]{\ref{lem:#1}}
\newcommand{\propref}[1]{Proposition~\ref{prop:#1}}
\newcommand{\propsref}[1]{Propositions~\ref{prop:#1}}
\newcommand{\propssref}[1]{\ref{prop:#1}}
\newcommand{\thmref}[1]{Theorem~\ref{thm:#1}}
\newcommand{\thmwref}[1]{Theorem~#1}
\newcommand{\thmsref}[1]{Theorems~\ref{thm:#1}}
\newcommand{\thmssref}[1]{\ref{thm:#1}}
\newcommand{\corref}[1]{Corollary~\ref{cor:#1}}
\newcommand{\tabref}[1]{Table~\ref{tab:#1}}
\newcommand{\tabsref}[1]{Tables~\ref{tab:#1}}
\newcommand{\tabssref}[1]{\ref{tab:#1}}
\newcommand{\pwref}[1]{p.~{#1}}
\newcommand{\ppwref}[1]{pp.~{#1}}
\newcommand{\eqnref}[1]{\eqref{eqn:#1}}




% \makeatletter
% \def\@strippedMR{}
% \def\@scanforMR#1#2#3\endscan{%
%   \ifx#1M\ifx#2R\def\@strippedMR{#3}%
%   \else\def\@strippedMR{#1#2#3}%
%   \fi\fi}
% \renewcommand\MR[1]{\relax
%   \ifhmode\unskip\spacefactor3000 \space\fi
%   \@scanforMR#1\endscan
%   MR\MRhref{\@strippedMR}{\@strippedMR}}
% \makeatother 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



\title[]{Order-invariant prior specification in Bayesian factor analysis}

\author[D.~Leung]{Dennis Leung} 
\address{Department of Statistics, University of Washington, Seattle,
  WA, U.S.A.}
\email{dmhleung@uw.edu}


\author[M.~Drton]{Mathias Drton} 
\address{Department of Statistics, University of Washington, Seattle,
  WA, U.S.A.}
\email{md5@uw.edu}



\begin{document}

\begin{abstract}
  In (exploratory) factor analysis, the loading matrix is identified
  only up to orthogonal rotation.  For identifiability, one thus often
  takes the loading matrix to be lower triangular with positive
  diagonal entries.  In Bayesian inference, a standard practice is
  then to specify a prior under which the loadings are independent,
  the off-diagonal loadings are normally distributed, and the diagonal
  loadings follow a truncated normal distribution.  This prior
  specification, however, depends in an important way on how the
  variables and associated rows of the loading matrix are ordered.  We
  show how a minor modification of the approach allows one to compute
  with the identifiable lower triangular loading matrix but maintain
  invariance properties under reordering of the variables.
\end{abstract}

\keywords{}

\subjclass[2000]{62H05}

\maketitle



\section{Introduction}
\label{sec:introduction}


Let $y$ be an $m$-vector of observed random variables, which for
simplicity we take to be centered.  Let $f\sim N_k(0,I_k)$ be a standard
normal $k$-vector of latent factors, with $k\le m$.  The factor
analysis model postulates that
\begin{equation} 
  \label{vecModel}
  {y} = {\beta}{f} + {\varepsilon}, 
\end{equation}
where $\beta=(\beta_{ij})\in\mathbb{R}^{m\times k}$ is an unknown
loading matrix, and $\varepsilon\sim N_m(0,\Omega)$ is an $m$-vector
of normally distributed error terms that are independent of $f$.  The
error terms are assumed to be mutually independent with
$\Omega=\text{diag}(\omega^2_1, \dots, \omega^2_m)$ comprising $m$
unknown positive variances that are also known as uniquenesses.  This
model with an unrestricted $m\times k$ loading matrix $\beta$ is
sometimes referred to as exploratory factor analysis---in contrast to
confirmatory factor analysis, which refers to situations in which some
collection of entries of $\beta$ is modeled as zero.

Integrating out the latent factors $f$ in~(\ref{vecModel}), the
observed random vector $y$ is seen to follow a centered multivariate
normal distribution with covariance matrix
\begin{equation}
  \label{eq:y:sigma}
  \Sigma = \Omega + \beta\beta'.
\end{equation}
As discussed in detail in \cite{MR0084943}, $\Sigma$ determines the
unrestricted loading matrix $\beta$ only up to orthogonal rotation.
Indeed, $\beta\beta'=\beta QQ'\beta'$ for any $k\times k$ orthogonal
matrix $Q$.  More details on factor analysis can be found, for
instance, in
\cite{MR2849614}, \cite{MR2299716},
and \cite{MR2604888}.

In this paper, we are concerned with Bayesian inference in
(exploratory) factor analysis.  In Bayesian computation, it is
convenient to impose an identifiability constraint on the loading
matrix $\beta$.  A common choice is to restrict $\beta$ to be lower
triangular with nonnegative diagonal entries, that is, $\beta_{ij}=0$
for $1\le i<j\le k$ and $\beta_{ii}\ge 0$ for $1\le i\le k$
\citep{geweke:1996,aguilar:2000,MR2036762}.  Under these constraints,
a full rank matrix $\beta$ is uniquely determined by $\beta\beta'$.
In the papers just referenced and also the software implementation
provided by \cite{mcmcpack}, a default prior on the lower triangular
loading matrix has all its non-zero entries independent with
\begin{align}
  \label{eq:beta:prior:standard}
  \beta_{ij} &\sim 
  \begin{cases}
    \mathit{TN}(0, C_0) &\text{ if } i=j,\\
    N(0, C_0) &\text{ if } i>j.
  \end{cases}
\end{align}
Here, $\mathit{TN}(0,C_0)$ denotes a truncated normal distribution on
$(0,\infty)$, i.e., the conditional distribution of $X$ given $X>0$
for $X\sim N(0,C_0)$.  The variance $C_0>0$ is a hyperparameter.  The
prior distribution for the uniquenesses has
$\omega_1^2,\dots,\omega_m^2$ independent of $\beta$ and also mutually
independent with Inverse Gamma distribution,
\begin{align}
  \label{eq:IGprior}
  \omega_i^2&\sim \mathit{IG}(\nu/2, \nu s^2/2)
\end{align}
for hyperparameters $\nu,s>0$.  Equivalently, $\nu s^2/\omega_i^2$ is
chi-square distributed with $\nu$ degrees of freedom; compare
Eqn.~(26) in
\cite{geweke:1996}.

As discussed in \citet[Sect.~6]{MR2036762}, the prior specification
in~(\ref{eq:beta:prior:standard}) is such that the induced prior on
$\beta\beta'$ and the covariance matrix $\Sigma$ in
(\ref{eq:y:sigma}) depends on the way the variables and the associated
rows of the loading matrix $\beta$ are ordered.  Indeed, a priori,
\begin{equation}
  \label{eq:betabeta:ii:chisquare}
  (\beta\beta')_{ii}/C_0=\sum_{j=1}^k
  \beta_{ij}^2/C_0=\sum_{j=1}^{\min\{i,k\}}
  \beta_{ij}^2/C_0
\end{equation}
follows a chi-square distribution with degrees of freedom
$\min\{i,k\}$.  Consequently, the implied prior and also the
posterior distribution for the covariance matrix $\Sigma$ is not
invariant under permutations of the variables.

In this paper we propose a modification of the prior distribution for
$\beta$ that maintains the convenience of computing with an
identifiable lower triangular loading matrix all the while making the
prior distributions of $\beta\beta'$ and $\Sigma$ invariant under
reordering of the variables.  Our proposal, described in
Section~\ref{sec:order-invariant}, merely changes the prior
distributions of the diagonal entries $\beta_{ii}$
in~(\ref{eq:beta:prior:standard}), which will be taken from a slightly
more general family than the truncated normal.  The details of a Gibbs
sampler to draw from the resulting posterior are given in
Section~\ref{sec:gibbs}.  We conclude with numerical examples and a
discussion in Sections~\ref{sec:numerics} and~\ref{sec:conclusion},
respectively.


\section{Order-invariant prior distribution}
\label{sec:order-invariant}

Without any identifiability constraints, the loading matrix $\beta$
takes its values in all of $\mathbb{R}^{m\times k}$.  A natural
default prior would then be to take all entries $\beta_{ij}$,
$i=1,\dots,m$ and $j=1,\dots, k$, to be independent $N(0,C_0)$ random
variables; we write $\beta\sim N_{m\times k}(0,C_0I_m\otimes I_k)$.
The spherical normal distribution $N_{m\times k}(0,C_0I_m\otimes I_k)$
is clearly invariant under permutation of the rows of the matrix.
Hence, the induced prior distribution of $\beta\beta'$ and of the
covariance matrix $\Sigma$ from~(\ref{eq:y:sigma}) is invariant under
simultaneous permutation of rows and columns.

Working with the prior just described comes at the cost of losing the
identifiability of $\beta$.  However, this can be overcome as follows.
Assuming that $m \geq k$, any $m \times k$ matrix $\beta$
with linearly independent columns can be uniquely decomposed as $\beta
= LQ$, where $L$ is an $m \times k$ lower triangular matrix with
positive diagonal, and $Q$ is a $k \times k$
orthogonal matrix.  We may then use the implied distribution of the
lower triangular matrix $L$ as a prior on the loading matrix.  The
following theorem about the joint distribution of $L$ and $Q$ is
adapted from Theorem $2.1.13$ in \cite{MR652932}.

\begin{theorem} 
  \label{LQthm} 
  Let $\beta=LQ$ be the LQ decomposition of the $m\times k$ random
  matrix $\beta\sim N_{m\times k}(0,C_0I_m\otimes I_k)$, where $m\ge k$.  Then the
  lower triangular matrix $L$ and the orthogonal matrix $Q$ are
  independent, the distribution of $Q$ is the normalized Haar measure,
  and the distribution of $L=(L_{ij})$ has joint density
  proportional to
  \begin{equation}
    \label{Lprior}
    \prod_{i=1}^m\prod_{j=1}^{\min\{i,k\}} \exp\left\{-\frac{1}{2C_0}L_{ij}^2 \right\}\times
    \prod_{i = 1}^k L_{ii}^{k - i} \mathbf{1}_{\{L_{ii}>0 \}}
  \end{equation}
  with respect to the Lebesgue measure on the space of $m \times k$
  lower triangular matrices.
\end{theorem}

The joint distribution for the entries of $L=(L_{ij})$ given
by~(\ref{Lprior}) has the entries $L_{ij}$, $i\ge j$, independent with
$L_{ij}\sim N(0,C_0)$ if $i>j$ and $L_{ii}$ following the distribution
with density proportional to
\begin{equation}
  \label{eq:new-distribution}
  x^{k-i}\exp\left\{-\frac{1}{2C_0}x^2\right\}, \quad x>0.
\end{equation}
Note that $L_{kk}\sim\mathit{TN}(0,C_0)$.  The joint distribution for
a lower triangular matrix in~(\ref{Lprior}) thus differs from that
given by~(\ref{eq:beta:prior:standard}) only in the coordinates
$L_{ii}$ for $1\le i\le k-1$, which are no longer truncated normal.

Assume as in~(\ref{eq:IGprior}) that $\Omega$ and $\beta$ are
independent a priori.  Then since $Q$ is independent of $L$, and
\[
\Sigma=\Omega+\beta\beta'=\Omega+LQQ'L' = \Omega+LL'
\]
does not depend on $Q$, the tuple $(y,\Omega,L,\Sigma)$ is independent
of $Q$.  Hence, $(\Omega,L,\Sigma)$ is also independent of $Q$ a
posteriori (i.e., conditional on $y$).  Our proposal is now simply to
keep with the standard identifiability constraint that has the loading
matrix $\beta$ lower triangular with nonnegative diagonal entries but to
use the distribution given by~(\ref{Lprior}) instead
of~(\ref{eq:beta:prior:standard}) for this lower triangular loading
matrix.  Concerning the remaining parts of the prior specification, we
continue to assume independence of $\beta$ and $\Omega$, and we stick
with the choice from~(\ref{eq:IGprior}) for the prior on the
uniquenesses.  This proposed prior has then the property that the
distributions of $\beta\beta'$ and the covariance matrix $\Sigma$ are
invariant under reordering of the variables (i.e., matrix rows and
columns), both a priori and a posteriori.

\section{Gibbs sampler}
\label{sec:gibbs}

Consider now an actual inferential setting in which we observe a
sample $y_1,\dots,y_n$ that comprises $n$ independent random vectors
drawn from a distribution in the $k$-factor model.  Let $Y$ be the
$n\times m$ matrix with the vectors $y_1,\dots,y_n$ as rows.  Let $F$
be an associated $n\times k$ matrix whose rows $f_1, \dots,f_n$ are
independent vectors of latent factors.  The factor analysis model
dictates that
\begin{equation} 
  \label{matModel}
  {Y} = {F} {\beta}' + E, 
\end{equation}
where $E= (\varepsilon_1, \dots, \varepsilon_n)'$ is an $n\times m$
matrix of stochastic errors.  The pairs $(f_t, \varepsilon_t)$ for
$1\le t\le n$ are independent, and in each pair $f_t\sim N_k(0,I_k)$
and $\varepsilon_t\sim N_m(0,\Omega)$ are independent as well.  The
unknown parameters are comprised in the matrices $\Omega =
\text{diag}(\omega^2_1, \dots, \omega^2_m)$ and $\beta =
(\beta_{ij})\in\mathbb{R}^{m\times k}$, where the latter is restricted to
be lower triangular with nonnegative diagonal.

We now adopt the prior distribution on $\beta$ and $\Omega$ given
by~(\ref{Lprior}) and~(\ref{eq:IGprior}), and derive the full
conditionals needed for a Gibbs sampler that generates draws from the
posterior distribution of $(\beta,\Omega)$.  As in \cite{MR2036762},
we write
\[
\beta_i =
\begin{cases}
  (\beta_{i1}, \dots, \beta_{ii})' &\text{ if } i\le k,\\
  (\beta_{i1}, \dots, \beta_{ik})' &\text{ if } i> k,
\end{cases}
\]
and explicitly involve the latent factors in $F$.  Let $F_i$ be the
$n\times i$ matrix made up of the first $i$ columns of $F$, and write
$Y_i$ for the $i$-th column of $Y$ (in contrast to $y_t$, which is the
$t$-th row of $Y$).  The full conditionals for $F$, $\Omega$ and
$\beta$ are determined as follows.  First, the rows $f_t$ of $F$ are
conditionally independent given $(\beta,\Omega,Y)$ with
\begin{equation}
  \label{eq:fullcond:F}
  (f_t \,|\,\beta, \Omega,Y) \;\sim\; N_k\left((I_k  + \beta'\Omega^{-1}\beta)^{-1} \beta' \Omega^{-1} y_t   , ({I}_k + {\beta}'{\Omega}^{-1}{\beta})^{-1}\right)
\end{equation}
for $t = 1, \dots, n$.  Second, the uniquenesses
$\omega_1^2,\dots,\omega_m^2$ are conditionally independent given
$(\beta,F,Y)$ with
\begin{equation}
  \label{eq:fullcond:omega}
(\omega^2_i \,|\,  \beta,F,Y) \;\sim\; \mathit{IG}\left(\tfrac{1}{2}(\nu + T),\tfrac{1}{2}(\nu
s^2 + d_i)\right),
\end{equation}
where 
\[
d_i = ({ Y}_i - {F}_i{\beta}_{i}' )'({Y}_i - {F}_i{\beta}_{i}').
\]
Third, the rows of $\beta$ are conditionally independent given $(\Omega,
F, Y)$.  For $i = 1, \dots, k$, the conditional density of the
vector $\beta_i$ is proportional to
\begin{equation}
  \label{eq:fullcond:beta:1}
   \beta_{ii}^{k - i}\frac{1}{\det(C_i)}\exp\left\{-\frac{1}{2}({\beta}_i - { m}_i)' {C}_i^{-1}({\beta}_i - {
     m}_i)\right\}\mathbf{ 1}_{\{\beta_{ii}>0\}},
\end{equation}
where 
\[
{C}_i
= \left(\frac{1}{C_0}{ I}_i + \frac{1}{\omega_i^{2}}{ F}_i'{
    F}_i\right)^{-1}
\quad\text{and}\quad 
{m}_i = \frac{1}{\omega_i^{2}}{C}_i   {F}_i' { Y}_i .
\]
For $i = k+1, \dots, m$, the conditional distribution is
\begin{equation}
  \label{eq:fullcond:beta:2}
  ({ \beta}_i\,|\,  \Omega,F,Y) \; \;\sim\; N_k({ m}_i, { C}_i)
\end{equation}
with
\[
{C}_i
= \left(\frac{1}{C_0}{ I}_k + \frac{1}{\omega_i^{2}}{ F}'{
    F}\right)^{-1}
\quad\text{and}\quad 
{m}_i = \frac{1}{\omega_i^{2}}{C}_i   {F}' { Y}_i .
\]

The only full conditional that differs from those given in
\cite{MR2036762} is the one for $\beta_i$ with $i\le k$
from~(\ref{eq:fullcond:beta:1}).  To draw from this distribution, we
first sample from $(\beta_{ii}\,|\, \Omega,F,Y)$ and then from
$(\beta_{i1},\dots,\beta_{ii-1}\,|\,\beta_{ii},\Omega,F,Y)$.  The
latter distribution is a multivariate normal distribution.  The only
new challenge is thus the sampling from $(\beta_{ii}\,|\,
\Omega,F,Y)$, which has density proportional to
\[
\beta_{ii}^{k - i} e^{-\frac{(\beta_{ii} - a)^2}{2 b^2}}{\bf
  1}_{\{\beta_{ii} > 0\}}
\]
for constants $a \in \mathbb{R}$ and $b>0$ determined by
$(\Omega,F,Y)$.  After scaling $\beta_{ii}$ by $b$, the problem
reduces to generating draws from distributions with density in the class
\begin{equation}
  \label{eq:new:distribution}
  f(x\,|\, \alpha,\gamma) = \frac{1}{Z(\alpha,\gamma)} x^{\alpha-1}
    e^{-(x-\gamma)^2}, \quad x>0,
\end{equation}
where $\alpha>0$ and $\gamma\in\mathbb{R}$ are two parameters, and
$Z(\alpha,\gamma)$ is the normalizing constant.  In the present
context, integer values of $\alpha$ are of interest.  The densities
in~(\ref{eq:new:distribution}) are log-concave, and we use adaptive
rejection sampling \citep{gilks1992adaptive} as implemented in the R
package \texttt{ars} to generate from them.


\section{Numerical experiments}
\label{sec:numerics}

We illustrate the use of the two different priors, obtained
from~(\ref{eq:beta:prior:standard}) and~(\ref{Lprior}), respectively,
on a simulated dataset $Y$ that involves $m=15$ variables and is of
size $n=30$.  The data are drawn from the $k=3$ factor distribution
given by the following loading matrix and uniquenesses:
 \newcommand{\sss}{\scriptstyle }
\begin{equation*}
  \beta_0 = \bordermatrix{
    & \sss 1 & \sss 2 & \sss 3 \cr
    \sss 1 & 0.97 &  0    & 0    \cr 
    \sss 2 &  0.04 &  0.90 & 0    \cr 
    \sss 3 &  1.00 & -1.12 & 0.57 \cr 
    \sss 4 &  2.03 &  0.42 & 0.57 \cr 
    \sss 5 &  0.31 &  0.47 & 0.09 \cr 
    \sss 6 &  0.43 & -0.21 & -0.35\cr 
    \sss 7 &  0.75 &  0.31 & 0.68 \cr 
    \sss 8 &  0.45 & -0.48 & -1.50\cr 
    \sss 9 & -2.21 &  1.45 & 0.38 \cr 
    \sss 10 &  1.98 & -0.30 & 0.96 \cr 
    \sss 11 & -2.63 &  0.41 & 1.09 \cr 
    \sss 12 & -0.72 &  1.39 & 0.97 \cr 
    \sss 13 & -0.88 &  2.01 & -0.39\cr 
    \sss 14 & -0.53 &  0.04 & 0.59 \cr 
    \sss 15 & -0.95 &  1.39 & 0.37 \cr 
  }, \qquad  
  \text{diag}(\Omega_0) = \bordermatrix{
    &\cr
    \sss 1 & 0.17 \cr 
    \sss 2 &  0.05 \cr 
    \sss 3 &  0.02 \cr 
    \sss 4 &  0.02 \cr 
    \sss 5 &  0.05 \cr 
    \sss 6 &  0.06 \cr 
    \sss 7 &  0.04 \cr 
    \sss 8 &  0.67 \cr 
    \sss 9 & -0.04 \cr 
    \sss 10 & 0.21 \cr 
    \sss 11 & 0.10 \cr 
    \sss 12 & 0.09 \cr 
    \sss 13 & 0.21 \cr 
    \sss 14 & 0.51 \cr 
    \sss 15 & 0.03 \cr 
  }.
\end{equation*}
We create a second data matrix $Y^\pi$ by permuting the columns of $Y$
based on the permutation $\pi$ from Table \ref{tablepi}, i.e. the
$i$-th column of $Y$ becomes the $\pi(i)$-th column of $Y^{\pi}$.  For
Bayesian inference, we choose the hyperparameters as \cite{MR2036762},
that is, $C_0 = 1$, $\nu= 2.2$ and $s = \sqrt{0.1/2.2}$.  Via Gibbs
sampling, we draw from the posterior distributions for
the covariance matrix $\Sigma=\Omega+\beta\beta'$ for each data set,
focusing on the factor analysis models $k=3$, and $k=6$ factors.  The
Gibbs samplers are initialized at the respective maximum likelihood
estimates for $(\beta, \Omega)$.  After a burnin of $10, 000$
iterations, we ran each sampler for $300, 000$ iterations.


\begin{table}[t]
\centering
\caption{The permutation $\pi$ used to reorder simulated data.}
\label{tablepi}
\begin{tabular}{c|rrrrrrrrrrrrrrr}
  \hline
  \hline
$i$ &   1 &   2 &   3 &   4 &   5 &   6 &   7 &   8 &   9 &  10 &  11 &  12 &  13 &  14 &  15 \\ 
\hline
  $\pi(i)$ &  10 &   14 &  13 &  15 &  12 &   6 &   7 &  2 &  11 &   9 &   8 &   3 &   5 &   1 &   4 \\ 
   \hline
\end{tabular}
\end{table}


Figures~\ref{figurek3} and~\ref{figurek6} show kernel density
estimates of the posterior densities of selected variances.  More
precisely, we compare the densities of $(\sigma_{ii}\,|\,Y)$ and
$(\sigma_{\pi(i),\pi(i)}\,|\,Y^\pi)$ for $i=1,8,14$.  Under our
proposed prior from~(\ref{Lprior}), the two posterior densities are
the same.  Indeed, the plots in the right hand columns of
Figures~\ref{figurek3} and~\ref{figurek6} show only minor
discrepancies due to Monte Carlo error.  The `standard prior'
from~(\ref{eq:beta:prior:standard}), however, results in visible
differences that are more pronounced for $k=6$, which is not
surprising as larger differences are possible among the degrees of
freedom of the chi-square prior for $(\beta\beta')_{ii}/C_0$;
recall~(\ref{eq:betabeta:ii:chisquare}).  Note that the observed
shifts in the posterior distributions under the `standard prior' are
explained by the different chi-square degrees of freedom.

\begin{figure}[t]
\centering
      \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k3_i1data_Self_prior_old}  
      \end{minipage} 
      \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k3_i1data_Self_prior_new}  
      \end{minipage}   
     
     \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k3_i8data_Self_prior_old}  
      \end{minipage}
 \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k3_i8data_Self_prior_new}  
      \end{minipage}

      \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k3_i14data_Self_prior_old}  
      \end{minipage}   
      \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k3_i14data_Self_prior_new}  
      \end{minipage}
      \caption{Posterior densities of $(\sigma_{ii}\,|\,Y)$, in black,
        and of $(\sigma_{\pi(i),\pi(i)}\,|\,Y^\pi)$, in grey, in
        factor analysis with $k=3$ factors, for $i=1,8,14$.
        The left column concerns the prior
        from~(\ref{eq:beta:prior:standard}), and the right column is
        based on the prior proposed in~(\ref{Lprior}).  }
\label{figurek3}
\end{figure}


\begin{figure}[t]
\centering
      \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k6_i1data_Self_prior_old}  
      \end{minipage} 
      \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k6_i1data_Self_prior_new}  
      \end{minipage}   
     
     \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k6_i8data_Self_prior_old}  
      \end{minipage}
 \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k6_i8data_Self_prior_new}  
      \end{minipage}

      \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k6_i14data_Self_prior_old}  
      \end{minipage}   
      \begin{minipage}{.48 \linewidth}
              \includegraphics[width=\linewidth]{density_postvar_k6_i14data_Self_prior_new}  
      \end{minipage}
\caption{Posterior densities of $(\sigma_{ii}\,|\,Y)$, in black, and
  of $(\sigma_{\pi(i),\pi(i)}\,|\,Y^\pi)$, in grey, in factor analysis
  with $k=6$ factors, for $i=1,8,14$.
  The left column concerns the
  prior from~(\ref{eq:beta:prior:standard}), and the right column is
  based on the prior proposed in~(\ref{Lprior}).  }
\label{figurek6}
\end{figure}






\section{Conclusion}
\label{sec:conclusion}

This paper proposes a prior distribution for the loading matrix in
factor analysis.  The proposal allows for computation with an
identifiable lower triangular loading matrix $\beta$ all the while
having the associated covariance matrix invariant under permutation of
the variables at hand.  The prior is intended as a possible default
when there is no reason to impose dependence among loadings or to
treat the loadings of different variables differently.  Concerning
possible departures from our default scenario, we remark that the
software of \cite{mcmcpack} also allows one to impose patterns of
zeros in the loading matrix $\beta$.  As mentioned earlier, the latter
situation is sometimes termed confirmatory factor analysis.  The
identifiability issues we addressed need not arise in that case as
orthogonal transformations will generally not preserve prescribed
zeros in the loading matrix.

Sampling from the posterior distribution resulting from the prior we
proposed is largely the same as for the `standard prior' that has been
used by several authors including 
\cite{geweke:1996} and \cite{MR2036762}.  The key difference is the
need to sample from distributions in the class specified
by~(\ref{eq:new:distribution}).  These distributions also appear in
the realm of multivariate $t$-distributions
\citep{finegold:aoas,finegold:ba}, although a square-root
transformation is necessary to match the setup there.  It thus seems
worthwhile to develop an efficient sampler targeting precisely this
family of distributions, which is a problem we are working on.

Finally, we emphasize that our proposal rests 
in an important way on
the fact that we derived it from a spherical joint normal distribution
for the loading matrix, namely, $\beta\sim N_{m\times
  k}(0,C_0I_m\otimes I_k)$.  Departures from this situation, even
merely including a non-zero mean for this matrix normal distribution,
seem to lead to a considerably more difficult scenario.


  \section*{Acknowledgments}

  This work was supported by the U.S.~National Science Foundation
  (DMS-1305154) and by the University of Washington Royalty Research
  Fund.


\bibliographystyle{ba}
\bibliography{fact_anal_bib}


\end{document}





factanal_v5.bbl
\begin{thebibliography}{12}
\newcommand{\enquote}[1]{``#1''}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{{\tt #1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi

\bibitem[{Aguilar and West(2000)}]{aguilar:2000}
Aguilar, O. and West, M. (2000).
\newblock \enquote{Bayesian dynamic factor models and variance matrix
  discounting for portfolio allocation.}
\newblock {\em Journal of Business \& Economic Statistics\/}, 18(3): 338--357.

\bibitem[{Anderson and Rubin(1956)}]{MR0084943}
Anderson, T.~W. and Rubin, H. (1956).
\newblock \enquote{Statistical inference in factor analysis.}
\newblock In {\em Proceedings of the {T}hird {B}erkeley {S}ymposium on
  {M}athematical {S}tatistics and {P}robability, 1954--1955, vol. {V}\/},
  111--150. University of California Press, Berkeley and Los Angeles.

\bibitem[{Bartholomew et~al.(2011)Bartholomew, Knott, and Moustaki}]{MR2849614}
Bartholomew, D., Knott, M., and Moustaki, I. (2011).
\newblock {\em Latent variable models and factor analysis\/}.
\newblock Wiley Series in Probability and Statistics. John Wiley \& Sons, Ltd.,
  Chichester, third edition.
\newblock A unified approach.

\bibitem[{Drton et~al.(2007)Drton, Sturmfels, and Sullivant}]{MR2299716}
Drton, M., Sturmfels, B., and Sullivant, S. (2007).
\newblock \enquote{Algebraic factor analysis: tetrads, pentads and beyond.}
\newblock {\em Probab. Theory Related Fields\/}, 138(3-4): 463--493.

\bibitem[{Finegold and Drton(2011)}]{finegold:aoas}
Finegold, M. and Drton, M. (2011).
\newblock \enquote{Robust graphical modeling with classical and alternative
  $t$-distributions.}
\newblock {\em Annals of Applied Statistics\/}, 5(2A): 1057--1080.

\bibitem[{Finegold and Drton(2014)}]{finegold:ba}
--- (2014).
\newblock \enquote{Robust {B}ayesian graphical modeling using {D}irichlet
  $t$-distributions.}
\newblock {\em Bayesian Analysis\/}, 9(3): 521--550.

\bibitem[{Geweke and Zhou(1996)}]{geweke:1996}
Geweke, J. and Zhou, G. (1996).
\newblock \enquote{Measuring the pricing error of the arbitrage pricing
  theory.}
\newblock {\em The Review of Financial Studies\/}, 9(2): 557--587.

\bibitem[{Gilks and Wild(1992)}]{gilks1992adaptive}
Gilks, W.~R. and Wild, P. (1992).
\newblock \enquote{Adaptive rejection sampling for {G}ibbs sampling.}
\newblock {\em Applied Statistics\/}, 41(2): 337--348.

\bibitem[{Lopes and West(2004)}]{MR2036762}
Lopes, H.~F. and West, M. (2004).
\newblock \enquote{Bayesian model assessment in factor analysis.}
\newblock {\em Statist. Sinica\/}, 14(1): 41--67.

\bibitem[{Martin et~al.(2011)Martin, Quinn, and Park}]{mcmcpack}
Martin, A.~D., Quinn, K.~M., and Park, J.~H. (2011).
\newblock \enquote{{MCMCpack}: {M}arkov Chain {M}onte {C}arlo in {R}.}
\newblock {\em Journal of Statistical Software\/}, 42(9): 22.

\bibitem[{Muirhead(1982)}]{MR652932}
Muirhead, R.~J. (1982).
\newblock {\em Aspects of multivariate statistical theory\/}.
\newblock John Wiley \& Sons, Inc., New York.
\newblock Wiley Series in Probability and Mathematical Statistics.

\bibitem[{Mulaik(2010)}]{MR2604888}
Mulaik, S.~A. (2010).
\newblock {\em Foundations of factor analysis\/}.
\newblock Statistics in the Social and Behavioral Sciences Series. CRC Press,
  Boca Raton, FL, second edition.

\end{thebibliography}
