#' Bayesian Factor Analysis
#'
#' @param formula a symbolic representation of the model to be
#' estimated, in the form \code{~ Y1 + Y2 + Y3}, where Y1, Y2, and Y3 are variables
#' of interest in factor analysis (manifest variables), assumed to be
#' normally distributed. The model requires a minimum of three manifest
#' variables contained in the
#' same dataset. The \code{+} symbol means ``inclusion'' not
#' ``addition.''
#' @param factors number of the factors to be fitted (defaults to 2).
#' @param model the name of a statistical model to estimate.
#' For a list of other supported models and their documentation see:
#' \url{http://docs.zeligproject.org/articles/}.
#' @param data the name of a data frame containing the variables
#' referenced in the formula or a list of multiply imputed data frames
#' each having the same variable names and row numbers (created by
#' \code{Amelia} or \code{\link{to_zelig_mi}}).
#' @param ... additional arguments passed to \code{zelig},
#' relevant for the model to be estimated.
#' @param by a factor variable contained in \code{data}. If supplied,
#' \code{zelig} will subset
#' the data frame based on the levels in the \code{by} variable, and
#' estimate a model for each subset. This can save a considerable amount of
#' effort. You may also use \code{by} to run models using MatchIt
#' subclasses.
#' @param cite If is set to 'TRUE' (default), the model citation will be printed
#' to the console.
#'
#' @details
#' In addition, \code{zelig()} accepts the following additional arguments for model specification:
#' \itemize{
#' \item \code{lambda.constraints}: list containing the equality or
#' inequality constraints on the factor loadings. Choose from one of the following forms:
#' \item \code{varname = list()}: by default, no constraints are imposed.
#' \item \code{varname = list(d, c)}: constrains the dth loading for the
#' variable named varname to be equal to c.
#' \item \code{varname = list(d, +)}: constrains the dth loading for the variable named varname to be positive;
#' \item \code{varname = list(d, -)}: constrains the dth loading for the variable named varname to be negative.
#' \item \code{std.var}: defaults to \code{FALSE} (manifest variables are rescaled to
#' zero mean, but retain observed variance). If \code{TRUE}, the manifest
#' variables are rescaled to be mean zero and unit variance.
#' }
#'
#' In addition, \code{zelig()} accepts the following additional inputs for \code{bayes.factor}:
#' \itemize{
#' \item \code{burnin}: number of the initial MCMC iterations to be discarded (defaults to 1,000).
#' \item \code{mcmc}: number of the MCMC iterations after burnin (defaults to 20,000).
#' \item \code{thin}: thinning interval for the Markov chain. Only every thin-th
#' draw from the Markov chain is kept. The value of mcmc must be divisible
#' by this value. The default value is 1.
#' \item \code{verbose}: defaults to FALSE. If TRUE, the
#' progress of the sampler (every 10%10%) is printed to the screen.
#' \item \code{seed}: seed for the random number generator. The default is NA which
#' corresponds to a random seed 12345.
#' \item \code{Lambda.start}: starting values of the factor loading matrix \eqn{\Lambda}, either a
#' scalar (all unconstrained loadings are set to that value), or a matrix with
#' compatible dimensions. The default is 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 \code{Psi.start}: starting values for the uniquenesses, either a scalar
#' (the starting values for all diagonal elements of \eqn{\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 \eqn{\Psi} take the values of
#' Psi.start. The default value is NA where the starting values of the all the
#' uniquenesses are set to be 0.5.
#' \item \code{store.lambda}: defaults to TRUE, which stores the posterior draws of the factor loadings.
#' \item \code{store.scores}: defaults to FALSE. If 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.)
#' }
#'
#' The model also accepts the following additional arguments to specify prior parameters:
#' \itemize{
#' \item \code{l0}: mean of the Normal prior for the factor loadings, either a scalar or a
#' matrix with the same dimensions as \eqn{\Lambda}. If a scalar value, that value will be the
#' prior mean for all the factor loadings. Defaults to 0.
#' \item \code{L0}: precision parameter of the Normal prior for the factor loadings, either
#' a scalar or a matrix with the same dimensions as \eqn{\Lambda}. If \code{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 \code{a0}: the shape parameter of the Inverse Gamma prior for the uniquenesses
#' is \code{a0}/2. It can take a scalar value or a vector. The default value is 0.001.
#' \item \code{b0}: the scale parameter of the Inverse Gamma prior for the uniquenesses
#' is \code{b0}/2. It can take a scalar value or a vector. The default value is 0.001.
#' }
#'
#' Additional parameters avaialable to this model include:
#' \itemize{
#' \item \code{weights}: vector of weight values or a name of a variable in the dataset
#' by which to weight the model. For more information see:
#' \url{http://docs.zeligproject.org/articles/weights.html}.
#' }
#'
#' @return Depending on the class of model selected, \code{zelig} will return
#' an object with elements including \code{coefficients}, \code{residuals},
#' and \code{formula} which may be summarized using
#' \code{summary(z.out)} or individually extracted using, for example,
#' \code{coef(z.out)}. See
#' \url{http://docs.zeligproject.org/articles/getters.html} for a list of
#' functions to extract model components. You can also extract whole fitted
#' model objects using \code{\link{from_zelig_model}}.
#' @examples
#' \dontrun{
#' data(swiss)
#' names(swiss) <- c("Fert", "Agr", "Exam", "Educ", "Cath", "InfMort")
#' z.out <- zelig(~ Agr + Exam + Educ + Cath + InfMort,
#' model = "factor.bayes", data = swiss,
#' factors = 2, verbose = FALSE,
#' a0 = 1, b0 = 0.15, burnin = 500, mcmc = 5000)
#'
#' z.out$geweke.diag()
#' z.out <- zelig(~ Agr + Exam + Educ + Cath + InfMort,
#' 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 = FALSE, a0 = 1, b0 = 0.15,
#' burnin = 500, mcmc = 5000)
#'
#' z.out$geweke.diag()
#' z.out$heidel.diag()
#' z.out$raftery.diag()
#' summary(z.out)
#' }
#'
#' @seealso Vignette: \url{http://docs.zeligproject.org/articles/zelig_factorbayes.html}
#' @import methods
#' @export Zelig-factor-bayes
#' @exportClass Zelig-factor-bayes
#'
#' @include model-zelig.R
zfactorbayes <- setRefClass("Zelig-factor-bayes",
contains = c("Zelig"))
zfactorbayes$methods(
initialize = function() {
callSuper()
.self$name <- "factor-bayes"
.self$year <- 2013
.self$authors <- "Ben Goodrich, Ying Lu"
.self$packageauthors <- "Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park"
.self$description = "Bayesian Factor Analysis"
.self$fn <- quote(MCMCpack::MCMCfactanal)
# JSON from parent
.self$wrapper <- "factor.bayes"
}
)
zfactorbayes$methods(
zelig = function(formula,
factors = 2,
burnin = 1000, mcmc = 20000,
verbose = 0,
...,
data,
by = NULL,
bootstrap = FALSE) {
if(!identical(bootstrap,FALSE)){
stop("Error: The bootstrap is not available for Markov chain Monte Carlo (MCMC) models.")
}
.self$zelig.call <- match.call(expand.dots = TRUE)
.self$model.call <- .self$zelig.call
if (missing(verbose))
verbose <- round((mcmc + burnin) / 10)
if (factors < 2)
stop("Number of factors needs to be at least 2")
.self$model.call$verbose <- verbose
.self$model.call$x <- formula
.self$model.call$factors <- factors
callSuper(formula = formula, data = data,..., by = by, bootstrap = FALSE)
}
)
zfactorbayes$methods(
qi = function() {
return(NULL)
}
)
# The following diagnostics are also in Zelig-bayes, which unfortunately Zelig-factor-bayes does not currently inherit.
zfactorbayes$methods(
geweke.diag = function() {
diag <- lapply(.self$zelig.out$z.out, coda::geweke.diag)
# Collapse if only one list element for prettier printing
if(length(diag)==1){
diag<-diag[[1]]
}
if(!citation("coda") %in% .self$refs){
.self$refs<-c(.self$refs,citation("coda"))
}
ref1<-bibentry(
bibtype="InCollection",
title = "Evaluating the accuracy of sampling-based approaches to calculating posterior moments.",
booktitle = "Bayesian Statistics 4",
author = person("John", "Geweke"),
year = 1992,
publisher = "Clarendon Press",
address = "Oxford, UK",
editor = c(person("JM", "Bernado"), person("JO", "Berger"), person("AP", "Dawid"), person("AFM", "Smith"))
)
.self$refs<-c(.self$refs,ref1)
return(diag)
}
)
zfactorbayes$methods(
heidel.diag = function() {
diag <- lapply(.self$zelig.out$z.out, coda::heidel.diag)
# Collapse if only one list element for prettier printing
if(length(diag)==1){
diag<-diag[[1]]
}
if(!citation("coda") %in% .self$refs){
.self$refs<-c(.self$refs,citation("coda"))
}
ref1<-bibentry(
bibtype="Article",
title = "Simulation run length control in the presence of an initial transient.",
author = c(person("P", "Heidelberger"), person("PD", "Welch")),
journal = "Operations Research",
volume = 31,
year = 1983,
pages = "1109--44")
.self$refs<-c(.self$refs,ref1)
return(diag)
}
)
zfactorbayes$methods(
raftery.diag = function() {
diag <- lapply(.self$zelig.out$z.out, coda::raftery.diag)
# Collapse if only one list element for prettier printing
if(length(diag)==1){
diag<-diag[[1]]
}
if(!citation("coda") %in% .self$refs){
.self$refs<-c(.self$refs,citation("coda"))
}
ref1<-bibentry(
bibtype="Article",
title = "One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo.",
author = c(person("Adrian E", "Raftery"), person("Steven M", "Lewis")),
journal = "Statistical Science",
volume = 31,
year = 1992,
pages = "1109--44")
ref2<-bibentry(
bibtype="InCollection",
title = "The number of iterations, convergence diagnostics and generic Metropolis algorithms.",
booktitle = "Practical Markov Chain Monte Carlo",
author = c(person("Adrian E", "Raftery"), person("Steven M", "Lewis")),
year = 1995,
publisher = "Chapman and Hall",
address = "London, UK",
editor = c(person("WR", "Gilks"), person("DJ", "Spiegelhalter"), person("S", "Richardson"))
)
.self$refs<-c(.self$refs,ref1,ref2)
return(diag)
}
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.