\SweaveOpts{results=hide, prefix.string=vigpics/factorBayes}
\include{zinput}
%\VignetteIndexEntry{Bayesian Factor Analysis}
%\VignetteDepends{Zelig, MCMCpack}
%\VignetteKeyWords{model,latent factors,continuous,explanatory variables, normal, Gibbs}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@

<<loadLibrary, echo=F,results=hide>>=
pkg <- search()
if(!length(grep("package:Zelig",pkg)))
library(Zelig)
@ 

\section{\texttt{factor.bayes}: Bayesian Factor Analysis}
\label{factor.bayes}

Given some unobserved explanatory variables and observed dependent
variables, the Normal theory factor analysis model estimates the
latent factors.  The model is implemented using a Markov Chain Monte
Carlo algorithm (Gibbs sampling with data augmentation).  For factor
analysis with ordinal dependent variables, see ordered factor analysis
(\Sref{factor.ord}), and for a mix of types of dependent variables,
see the mixed factor analysis model (\Sref{factor.mix}).

\subsubsection{Syntax}
\begin{verbatim}
> z.out <- zelig(cbind(Y1 ,Y2, Y3) ~ NULL, factors = 2, 
                model = "factor.bayes", data = mydata)
\end{verbatim}

\subsubsection{Inputs}
{\tt zelig()} takes the following functions for {\tt factor.bayes}:  
\begin{itemize}
\item \texttt{Y1}, {\tt Y2}, and \texttt{Y3}: variables of interest in 
factor analysis (manifest variables), assumed to be normally
distributed. The model requires a minimum of three manifest variables.

\item \texttt{factors}: number of the factors to be fitted (defaults to 2).

\end{itemize}

\subsubsection{Additional Inputs}

In addition, {\tt zelig()} accepts the following additional arguments for 
model specification: 
\begin{itemize}
\item \texttt{lambda.constraints}: list containing the equality or 
inequality constraints on the factor loadings.  Choose from one of the
following forms:
\begin{itemize}

\item {\tt varname = list()}: by default, no constraints are
imposed.

\item \texttt{varname = list(d, c)}: constrains the
$d$th loading for the variable named \texttt{varname} to be equal to \texttt{c}.

\item \texttt{varname = list(d, "+")}: constrains the
$d$th loading for the variable named \texttt{varname} to be positive;

\item \texttt{varname = list(d, "-")}: constrains the
$d$th loading for the variable named \texttt{varname} to be negative.
\end{itemize} 

\item \texttt{std.var}: defaults to {\tt FALSE} (manifest variables
are rescaled to zero mean, but retain observed variance).  If
\texttt{TRUE}, the manifest variables are rescaled to be mean zero and
unit variance.

\end{itemize}
 In addition, {\tt zelig()} accepts the following additional inputs
for {\tt bayes.factor}:

\begin{itemize}
\item \texttt{burnin}: number of the initial MCMC iterations to be 
 discarded (defaults to 1,000). 

\item \texttt{mcmc}: number of the MCMC iterations after burnin
(defaults to 20,000).

\item \texttt{thin}: thinning interval for the Markov chain. Only every 
 \texttt{thin}-th draw from the Markov chain is kept. The value of 
\texttt{mcmc} must be divisible by this value. The default value is 1.

\item \texttt{verbose}: defaults to {\tt FALSE}.  If \texttt{TRUE}, the progress 
 of the sampler (every $10\%$) is printed to the screen.

\item \texttt{seed}: seed for the random number generator. The default
is \texttt{NA} which corresponds to a random seed 12345. 

\item \texttt{Lambda.start}: starting values of the factor loading
matrix $\Lambda$, either a scalar (all unconstrained loadings are set
to that value), or a matrix with compatible dimensions.  The default
is \texttt{NA}, where the start value are set to be 0 for
unconstrained factor loadings, and 0.5 or $-$0.5 for constrained
factor loadings (depending on the nature of the constraints).

\item \texttt{Psi.start}: starting values for the uniquenesses, either
a scalar (the starting values for all diagonal elements of $\Psi$ are
set to be this value), or a vector with length equal to the number of
manifest variables.  In the latter case, the starting values of the
diagonal elements of $\Psi$ take the values of \texttt{Psi.start}. The
default value is \texttt{NA} where the starting values of the all the
uniquenesses are set to be 0.5.

\item \texttt{store.lambda}: defaults to {\tt TRUE}, which stores the
posterior draws of the factor loadings.  

\item \texttt{store.scores}: defaults to {\tt FALSE}.  If {\tt TRUE},
stores the posterior draws of the factor scores.  (Storing factor
scores may take large amount of memory for a large number of draws
or observations.)

\end{itemize}

\noindent The model also accepts the following additional arguments to 
specify prior parameters:

\begin{itemize}

\item \texttt{l0}: mean of the Normal prior for the factor
loadings, either a scalar or a matrix with the same dimensions as
$\Lambda$.  If a scalar value, that value will be the prior mean for
all the factor loadings. Defaults to 0.

\item \texttt{L0}: precision parameter of the Normal prior 
for the factor loadings, either a scalar or a matrix with the same
dimensions as $\Lambda$.  If \texttt{L0} takes a scalar value, then
the precision matrix will be a diagonal matrix with the diagonal
elements set to that value. The default value is 0, which leads to an
improper prior.

\item \texttt{a0}: the shape parameter of the Inverse Gamma prior for 
the uniquenesses is \texttt{a0/2}. It can take a scalar value or a
vector. The default value is 0.001.

\item \texttt{b0}: the shape parameter of the Inverse Gamma prior for 
the uniquenesses is \texttt{b0/2}. It can take a scalar value or a
vector. The default value is 0.001.

\end{itemize}

Zelig users may wish to refer to \texttt{help(MCMCfactanal)} for more 
information.

\input{coda_diag}

\subsubsection{Examples}

\begin{enumerate}
\item {Basic Example} \\
Attaching the sample  dataset:
<<BasicExample.data>>=
 data(swiss)
 names(swiss) <- c("Fert","Agr","Exam","Educ","Cath","InfMort")
@ 
Factor analysis:
<<BasicExample.zelig>>=
 z.out <- zelig(cbind(Agr, Exam, Educ, Cath, InfMort) ~ NULL, 
                 model = "factor.bayes", data = swiss, factors = 2,
	         verbose = TRUE, a0 = 1, b0 = 0.15, 
                 burnin = 5000, mcmc = 50000)
@ 
%it is not in the demo and it fails 
Checking for convergence before summarizing the estimates:
<<BasicExample.coeff>>=
 algor <- try(geweke.diag(z.out$coefficients), silent=T)
if(class(algor)=="try-error")
print(algor)
@  

Since the algorithm did not converge, we now add some constraints on
$\Lambda$. 

\item {Putting Constraints on $\Lambda$} \\
Put constraints on Lambda to optimize the algorithm:
<<Constraints.zelig>>= 
z.out <- zelig(cbind(Agr, Exam, Educ, Cath, InfMort) ~ NULL,  
               model = "factor.bayes", data = swiss, factors = 2,
               lambda.constraints = list(Exam=list(1,"+"),
                                    Exam=list(2,"-"), Educ=c(2,0),
                                    InfMort=c(1,0)), 
               verbose = TRUE, a0 = 1, b0 = 0.15, 
               burnin = 5000, mcmc = 50000)

 geweke.diag(z.out$coefficients)
 heidel.diag(z.out$coefficients)
 raftery.diag(z.out$coefficients)
 summary(z.out)
@  
\end{enumerate}

\subsubsection{Model} 
Suppose for observation $i$ we observe $K$ variables and hypothesize
that there are $d$ underlying factors such that:  
\begin{eqnarray*}
Y_i = \Lambda \phi_i+\epsilon_i
\end{eqnarray*}
where $Y_{i}$ is the vector of $K$ manifest variables for observation
$i$. $\Lambda$ is the $K \times d$ factor loading matrix and $\phi_i$
is the $d$-vector of latent factor scores. Both $\Lambda$ and $\phi$
need to be estimated.

\begin{itemize}
\item The \emph{stochastic component} is given by:
\begin{eqnarray*}
\epsilon_{i}  \sim \textrm{Normal}(0, \Psi).
\end{eqnarray*}
where $\Psi$ is a diagonal, positive definite matrix. The diagonal elements
of $\Psi$ are referred to as uniquenesses.

\item The \emph{systematic component} is given by
\begin{eqnarray*}
\mu_i = E(Y_i) = \Lambda\phi_i
\end{eqnarray*}

\item The independent conjugate \emph{prior} for each $\Lambda_{ij}$  is given by
\begin{eqnarray*}
\Lambda_{ij} \sim \textrm{Normal}(l_{0_{ij}}, L_{0_{ij}}^{-1})
\textrm{ for } i=1,\ldots, k; \quad j=1,\ldots, d. 
\end{eqnarray*}

\item The independent conjugate \emph{prior} for each $\Psi_{ii}$ is given by
\begin{eqnarray*}
\Psi_{ii} \sim \textrm{InverseGamma}(\frac{a_0}{2}, \frac{b_0}{2}), \textrm{ for } 
i = 1, \ldots, k.
\end{eqnarray*}

\item The \emph{prior} for $\phi_i$ is
\begin{eqnarray*}
\phi_i &\sim& \textrm{Normal}(0, I_d), \textrm{ for } i = 1, \ldots, n.
\end{eqnarray*}
where $I_d$ is a $ d\times d $ identity matrix.
\end{itemize}

\subsubsection{Output Values}

The output of each Zelig command contains useful information which you may
view. For example, if you run:
\begin{verbatim}
z.out <- zelig(cbind(Y1, Y2, Y3), model = "factor.bayes", data)
\end{verbatim}

\noindent then you may examine the available information in \texttt{z.out} by
using \texttt{names(z.out)}, see the draws from the posterior distribution of
the \texttt{coefficients} by using \texttt{z.out\$coefficients}, and view a default
summary of information through \texttt{summary(z.out)}. Other elements
available through the \texttt{\$} operator are listed below.

\begin{itemize}
\item From the \texttt{zelig()} output object \texttt{z.out}, you may extract:

\begin{itemize}
\item \texttt{coefficients}: draws from the posterior distributions
of the estimated factor loadings and the uniquenesses.  If
\texttt{store.scores = TRUE}, the estimated factors scores are also
contained in \texttt{coefficients}.

\item \texttt{data}: the name of the input data frame.
\item \texttt{seed}: the random seed used in the model.   

\end{itemize}

\item Since there are no explanatory variables, the \texttt{sim()}
procedure is not applicable for factor analysis models.

\end{itemize}


\subsection* {How to Cite} 

\input{cites/factor.bayes}
\input{citeZelig}

<<afterpkgs, echo=FALSE>>=
 after<-search()
 torm<-setdiff(after,before)
 for (pkg in torm)
 detach(pos=match(pkg,search()))
@

 \end{document}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:

