# R/faMain.R In fungible: Psychometric Functions from the Waller Lab

#### Documented in faMain

#' Automatic Factor Rotation from Random Configurations with Bootstrap Standard Errors
#'
#' This function conducts factor rotations (using the \pkg{GPArotation} package)
#' from a user-specified number of random (orthogonal) starting configurations.
#' Based on the resulting complexity function value, the function determines the
#' number of local minima and, among these local solutions, will find the
#' "global minimum" (i.e., the minimized complexity value from the finite
#' number of solutions). See Details below for an elaboration on the global
#' minimum. This function can also return bootstrap standard errors of the factor solution.
#'
#' @param X (Matrix) A raw data matrix (or data frame).
#' @param R (Matrix) A correlation matrix.
#' @param n (Numeric) Sample size associated with the correlation matrix. Defaults to n = NULL.
#' @param numFactors (Numeric) The number of factors to extract for subsequent rotation.
#' @param facMethod (Character) The method used for factor extraction
#' (\code{\link{faX}}). The supported options are "fals" for unweighted least
#' squares, "faml" for maximum likelihood, "fapa" for iterated principal axis
#' factoring, "faregLS" for regularized least squares,
#' "faregML" for regularized maximum likelihood, and "pca" for principal components
#'  analysis. The default method  is "fals".
#' \itemize{
#'   \item \strong{"fals"}: Factors are extracted using the unweighted least
#'   squares estimation procedure using the \code{\link{fals}} function.
#'   \item \strong{"faml"}: Factors are extracted using the maximum likelihood
#'   estimation procedure using the \code{\link[stats]{factanal}} function.
#'   \item \strong{"fapa"}: Factors are extracted using the iterated principal
#'   axis factoring estimation procedure using the \code{\link{fapa}} function.
#'   \item \strong{"faregLS"}: Factors are extracted using regularized
#'   least squares factor analysis using the \code{\link{fareg}} function.
#'   \item \strong{"faregML"}: Factors are extracted using regularized
#'   maximum likelihood factor using the \code{\link{fareg}} function.
#'   \item \strong{"pca"}: Principal components are extracted.
#' }
#' @param rotate (Character) Designate which rotation algorithm to apply. The
#' following are available rotation options: "oblimin", "quartimin", "targetT",
#' "targetQ", "oblimax", "entropy", "quartimax", "varimax", "simplimax",
#' "bentlerT", "bentlerQ", "tandemI", "tandemII", "geominT", "geominQ", "cfT",
#' "cfQ", "infomaxT", "infomaxQ", "mccammon", "bifactorT", "bifactorQ", and
#' "none". Defaults to rotate = "oblimin". See \pkg{GPArotation} package for more
#' details. Note that rotations ending in "T" and "Q" represent orthogonal and
#' oblique rotations, respectively.
#' @param targetMatrix (Matrix) This argument serves two functions. First, if a
#' user has requested either a "targetT" or "targetQ' rotation, then
#'  the target matrix is used to conduct a fully or partially
#' specified target rotation. In the latter case,  freely estimated factor
#' loadings are designated by "NA" values and rotation will be conducted using
#' Browne's (1972a, 1972b, 2001) method for a partially-specified
#' target rotation. Second, if any other rotation option is chosen then all
#' (but not rotated) with the target solution.
#' @param bootstrapSE (Logical) Computes bootstrap standard errors. All bootstrap
#' samples are aligned to the global minimum solution. Defaults to
#' bootstrapSE = FALSE (no standard errors).
#' @param numBoot (Numeric) The number bootstraps. Defaults to numBoot = 1000.
#' @param CILevel (Numeric) The confidence level (between 0 and 1) of the bootstrap
#' confidence interval. Defaults to CILevel = .95.
#' @param Seed (Numeric) Starting seed for reproducible bootstrap results and factor rotations.
#' Defaults to Seed = 1.
#' @param digits (Numeric) Rounds the values to the specified number of decimal
#' places. Defaults to digits = NULL (no rounding).
#' @param faControl (List) A list of optional parameters passed to the factor
#' \itemize{
#'   \item \strong{treatHeywood}: (Logical) In \code{fals}, if treatHeywood is
#'   true, a penalized least squares function is used to bound the communality
#'   estimates below 1.0. Defaults to treatHeywood = TRUE.
#'   \item \strong{nStart}: (Numeric) The number of starting values to be tried
#'   in \code{faml}. Defaults to nStart = 10.
#'   \item \strong{start}: (Matrix) NULL or a matrix of starting values, each column
#'   giving an initial set of uniquenesses. Defaults to start = NULL.
#'   \item \strong{maxCommunality}: (Numeric) In \code{faml}, set the maximum
#'   communality value for the estimated solution. Defaults to maxCommunality = .995.
#'   \item \strong{epsilon}: (Numeric) In \code{fapa}, the numeric threshold
#'   designating when the algorithm has converged. Defaults to epsilon = 1e-4.
#'   \item \strong{communality}: (Character) The method used to estimate the
#'   initial communality values in \code{fapa}. Defaults to communality = 'SMC'.
#'   \itemize{
#'     \item \strong{"SMC"}: Initial communalities are estimated by taking the
#'     squared multiple correlations of each indicator after regressing the
#'     indicator on the remaining variables.
#'     \item \strong{"maxr"}: Initial communalities equal the largest
#'     (absolute value) correlation in each column of the correlation matrix.
#'     \item \strong{"unity"}: Initial communalities equal 1.0 for all variables.
#'   }
#'   \item \strong{maxItr}: (Numeric) In \code{fapa}, the maximum number of
#'   iterations to reach convergence. Defaults to maxItr = 15,000.
#' }
#'
#' @param rotateControl (List) A list of control values to pass to the factor rotation algorithms.
#' \itemize{
#'   \item \strong{numberStarts}: (Numeric) The number of random (orthogonal)
#'   starting configurations for the chosen rotation method (e.g., oblimin). The first
#'   rotation will always commence from the unrotated factors orientation.
#'   Defaults to numberStarts = 10.
#'   \item \strong{gamma}: (Numeric) This is a tuning parameter (between 0
#'   and 1, inclusive) for an oblimin rotation.  See the \pkg{GPArotation}
#'   library's oblimin documentation for more details. Defaults to gamma = 0
#'   (i.e., a quartimin rotation).
#'   \item \strong{delta}: (Numeric) This is a tuning parameter for the geomin
#'    rotation. It adds a small number (default = .01) to the squared factor
#'   \item \strong{kappa}: (Numeric) The main parameterization of the
#'   Crawford-Ferguson (CF) rotations (i.e., "cfT" and "cfQ" for orthogonal and
#'   oblique CF rotation, respectively). Defaults to kappa = 0.
#'   \item \strong{k}: (Numeric) A specific parameter of the simplimax rotation.
#'   Defaults to k = the number of observed variables.
#'   \item \strong{standardize}: (Character) The standardization routine used
#'   on the unrotated factor structure. The three options are "none", "Kaiser",
#'   and "CM". Defaults to standardize = "none".
#'   \itemize{
#'     \item \strong{"none"}: No standardization is applied to the unrotated
#'     factor structure.
#'     \item \strong{"Kaiser"}: Use a factor structure matrix that has been
#'     normed by Kaiser's method (i.e., normalize all rows to have a unit length).
#'     \item \strong{"CM"}: Use a factor structure matrix that has been normed
#'      by the Cureton-Mulaik method.
#'   }
#'   \item \strong{epsilon}: (Numeric) The rotational convergence criterion to
#'   use. Defaults to epsilon = 1e-5.
#'   in the \code{\link{promaxQ}} rotation. Defaults to power = 4.
#'   \item \strong{maxItr}: (Numeric) The maximum number of iterations for the
#'   rotation algorithm. Defaults to maxItr = 15000.
#' }
#'
#' @param ... Values to be passed to the \code{\link[stats]{cor}} function.
#' \itemize{
#'   \item \strong{use}: (Character) A character string giving a method for
#'   computing correlations in the presence of missing values: "everything"
#'   (the default), "all.obs", "complete.obs", "na.or.complete", or
#'   "pairwise.complete.obs".
#'   \item \strong{method}: (Character) A character string indicating which
#'   correlation coefficient is to be computed: "pearson" (the default),
#'   "kendall", or "spearman".
#'   \item \strong{na.rm}: (Logical) Should missing values be removed (TRUE)
#'   or not (FALSE)?
#' }
#'
#' @details
#' \itemize{
#'   \item \strong{Global Minimum}: This function uses several random starting
#'   configurations for factor rotations in an attempt to find the global
#'   minimum solution. However, this function is not guaranteed to find the
#'   global minimum. Furthermore, the global minimum solution need not be
#'   more psychologically interpretable than any of the local solutions (cf.
#'   Rozeboom, 1992). As is recommended, our function returns all local
#'   solutions so users can make their own judgements.
#'   \item \strong{Finding clusters of local minima}: We find local-solution sets by sorting the rounded
#'   rotation complexity values (to the number of  digits specified in the \code{epsilon}
#'   argument of the \code{rotateControl} list) into sets with equivalent values. For example,
#'   by default \code{epsilon = 1e-5.} and thus \code{} will only evaluate the complexity
#'   values to five significant digits. Any differences beyond that value will not effect the final sorting.
#' }
#'
#' @return The \code{faMain} function will produce a lot of output in addition
#' to the rotated factor pattern matrix and the factor correlations.
#' \itemize{
#'   \item \strong{R}: (Matrix) Returns the correlation matrix, useful when raw data are supplied.
#'   lowest evaluated discrepancy function. This solution has the lowest
#'   discrepancy function \emph{of the examined random starting configurations}.
#'   It is not guaranteed to find the "true" global minimum. Note that multiple
#'    (or even all) local solutions can have the same discrepancy functions.
#'   \item \strong{Phi}: (Matrix) The factor correlations of the rotated factor
#'   solution with the lowest evaluated discrepancy function (see Details).
#'   \item \strong{facIndeterminacy}: (Vector) A vector (with length equal to the number of factors)
#'   containing Guttman's (1955) index of factor indeterminacy for each factor.
#'   \item \strong{h2}: (Vector) The vector of final communality estimates.
#'   errors across the bootstrapped factor solutions. Each matrix element is
#'   the standard deviation of all bootstrapped factor loadings for that element position.
#'   \item \strong{CILevel} (Numeric) The user-defined confidence level (between 0 and 1) of the bootstrap
#'    confidence interval. Defaults to CILevel = .95.
#'   interval width is specified by the user.
#'   interval width is specified by the user.
#'   \item \strong{PhiSE}: (Matrix) The matrix of factor correlation standard
#'   errors across the bootstrapped factor solutions. Each matrix element is
#'   the standard deviation of all bootstrapped factor correlations for that element position.
#'   \item \strong{PhiCIupper}: (Matrix) Contains the upper confidence interval
#'   of the bootstrapped factor correlation matrix. The confidence interval
#'   width is specified by the user.
#'   \item \strong{PhiCIlower}: (Matrix) Contains the lower confidence interval
#'   of the bootstrapped factor correlation matrix. The confidence interval
#'   width is specified by the user.
#'   \item \strong{facIndeterminacySE}: (Matrix) A row vector containing the
#'   standard errors of Guttman's (1955) factor indeterminacy indices across the
#'   bootstrap factor solutions.
#'   \item \strong{localSolutions}: (List) A list containing all local solutions
#'   in ascending order of their factor loadings, rotation complexity values (i.e., the first solution
#'   is the "global" minimum). Each solution returns the
#'   \itemize{
#'      \item \strong{Phi}: (Matrix) factor correlations,
#'      \item \strong{RotationComplexityValue}: (Numeric) the complexity value of the rotation algorithm,
#'      \item \strong{facIndeterminacy}: (Vector) A vector of factor indeterminacy indices for each common factor, and
#'      \item \strong{RotationConverged}: (Logical) convergence status of the rotation algorithm.
#'      }
#'   \item \strong{numLocalSets} (Numeric) How many sets of local solutions
#'    with the same discrepancy value were obtained.
#'   \item \strong{localSolutionSets}: (List) A list containing the sets of
#'   unique local minima solutions. There is one list element for every unique
#'   local solution that includes (a) the factor loadings matrix, (b) the factor
#'   correlation matrix (if estimated), and (c) the discrepancy value of the rotation algorithm.
#'   number of bootstrapped samples (representing the row, column, and depth, respectively).
#'   \item \strong{PhiArray}: (Array) Contains an array of all bootstrapped
#'   factor correlations. The dimension are the number of factors, the number
#'   of factors, and the number of bootstrapped samples (representing the row,
#'    column, and depth, respectively).
#'   \item \strong{facIndeterminacyArray}: (Array) Contains an array of all
#'   bootstrap factor indeterminacy indices. The dimensions are 1, the number
#'   of factors, and the number of bootstrap samples (representing the row,
#'   column, and depth order, respectively).
#'   \item \strong{faControl}: (List) A list of the control parameters passed
#'   to the factor extraction (\code{\link{faX}}) function.
#'   \item \strong{faFit}: (List) A list of additional output from the factor
#'   extraction routines.
#'  \itemize{
#'     \item \strong{facMethod}: (Character) The factor extraction routine.
#'     \item \strong{df}: (Numeric) Degrees of Freedom from the maximum
#'     likelihood factor extraction routine.
#'     \item \strong{n}: (Numeric) Sample size associated with the correlation matrix.
#'     \item \strong{objectiveFunc}: (Numeric) The evaluated objective function for the
#'     maximum likelihood factor extraction routine.
#'     \item \strong{RMSEA}: (Numeric) Root mean squared error of approximation
#'     from Steiger & Lind (1980). Note that bias correction is computed if the
#'     sample size is provided.
#'     \item \strong{testStat}: (Numeric) The significance test statistic for the maximum
#'     likelihood procedure. Cannot be computed unless a sample size is provided.
#'     \item \strong{pValue}: (Numeric) The p value associated with the significance test
#'     statistic for the maximum likelihood procedure. Cannot be computed unless
#'     a sample size is provided.
#'     \item \strong{gradient}: (Matrix) The solution gradient for the least squares factor
#'     extraction routine.
#'     \item \strong{maxAbsGradient}: (Numeric) The maximum absolute value of the
#'     gradient at the least squares solution.
#'     \item \strong{Heywood}: (Logical) TRUE if a Heywood case was produced.
#'     \item \strong{convergedX}: (Logical) TRUE if the factor
#'     \strong{extraction} routine converged.
#'     \item \strong{convergedR}: (Logical) TRUE if the factor \strong{rotation}
#'     routine converged (for the local solution with the minimum discrepancy
#'     value).
#'   }
#'   \item \strong{rotateControl}: (List) A list of the control parameters
#'   passed to the rotation algorithm.
#'   \item \strong{unSpunSolution}: (List) A list of output parameters (e.g., loadings, Phi, etc) from
#'   the rotated solution that was obtained by rotating directly from the unrotated (i.e., unspun) common factor orientation.
#'   \item \strong{targetMatrix} (Matrix) The input target matrix if supplied by the user.
#'   \item \strong{Call}: (call) A copy of the function call.
#' }
#'
#' @references Browne, M. W.  (1972).  Oblique rotation to a partially specified
#' target.  \emph{British Journal of Mathematical and Statistical Psychology, 25},(1),
#' 207-212.
#' @references Browne, M. W. (1972b). Orthogonal rotation to a partially specifed
#' target. \emph{British Journal of Statistical Psychology, 25},(1), 115-120.
#' @references Browne, M. W. (2001). An overview of analytic rotation in
#' exploratory factor analysis. \emph{Multivariate Behavioral Research, 36}(1), 111-150.
#' @references Cureton, E. E., & Mulaik, S. A. (1975). The weighted varimax
#' rotation and the promax rotation. \emph{Psychometrika, 40}(2), 183-195.
#' @references Guttman, L. (1955). The determinacy of factor score matrices with
#' implications for five other basic problems of common factor theory.
#' \emph{British Journal of Statistical Psychology, 8}(2), 65-81.
#' @references Jung, S. & Takane, Y.  (2008).  Regularized common factor analysis.
#' \emph{New Trends in Psychometrics}, 141-149.
#' @references Mansolf, M., & Reise, S. P. (2016). Exploratory bifactor
#' analysis: The Schmid-Leiman orthogonalization and Jennrich-Bentler
#' analytic rotations. \emph{Multivariate Behavioral Research, 51}(5), 698-717.
#' @references Rozeboom, W. W. (1992). The glory of suboptimal factor rotation:
#' Why local minima in analytic optimization of simple structure are more
#' blessing than curse. \emph{Multivariate Behavioral Research, 27}(4), 585-599.
#' @references Zhang, G. (2014). Estimating standard errors in exploratory factor
#' analysis. \emph{Multivariate Behavioral Research, 49}(4), 339-353.
#'
#' @family Factor Analysis Routines
#'
#' @author
#' \itemize{
#'   \item Niels G. Waller (nwaller@umn.edu)
#'   \item Casey Giordano (Giord023@umn.edu)
#'   \item The authors thank Allie Cooperman and Hoang
#'    Nguyen for their help implementing the standard error estimation and the
#'    Cureton-Mulaik standardization procedure.
#'}
#'
#' @examples
#' ## Example 1
#'
#' ## Generate an orthogonal factor model
#' lambda <- matrix(c(.41, .00, .00,
#'                    .45, .00, .00,
#'                    .53, .00, .00,
#'                    .00, .66, .00,
#'                    .00, .38, .00,
#'                    .00, .66, .00,
#'                    .00, .00, .68,
#'                    .00, .00, .56,
#'                    .00, .00, .55),
#'                  nrow = 9, ncol = 3, byrow = TRUE)
#'
#' ## Generate factor correlation matrix
#' Phi <- matrix(.50, nrow = 3, ncol = 3)
#' diag(Phi) <- 1
#'
#' ## Model-implied correlation matrix
#' R <- lambda %*% Phi %*% t(lambda)
#' diag(R) <- 1
#'
#' ## Load the MASS package to create multivariate normal data
#' library(MASS)
#'
#' ## Generate raw data to perfectly reproduce R
#' X <- mvrnorm(Sigma = R, mu = rep(0, nrow(R)), empirical = TRUE, n = 300)
#'
#'\dontrun{
#' ## Execute 50 promax rotations from a least squares factor extraction
#' ## Compute 100 bootstrap samples to compute standard errors and
#' ## 80 percent confidence intervals
#' Out1 <- faMain(X             = X,
#'                numFactors    = 3,
#'                facMethod     = "fals",
#'                rotate        = "promaxQ",
#'                bootstrapSE   = TRUE,
#'                numBoot       = 100,
#'                CILevel       = .80,
#'                faControl     = list(treatHeywood = TRUE),
#'                rotateControl = list(numberStarts = 2,
#'                                     power        = 4,
#'                                     standardize  = "Kaiser"),
#'                digits        = 2)
#'}
#'
#' ## Example 2
#'
#' ## Load Thurstone's (in)famous box data
#' data(Thurstone, package = "GPArotation")
#'
#' ## Execute 5 oblimin rotations with Cureton-Mulaik standardization
#'                rotate        = "oblimin",
#'                bootstrapSE   = FALSE,
#'                rotateControl = list(numberStarts = 5,
#'                                     standardize  = "CM",
#'                                     gamma        = 0,
#'                                     epsilon      = 1e-6),
#'                digits        = 2)
#'
#'
#' ## Example 3
#'
#' ## Factor matrix from Browne 1972
#' lambda <- matrix(c(.664,  .322, -.075,
#'                    .688,  .248,  .192,
#'                    .492,  .304,  .224,
#'                    .837, -.291,  .037,
#'                    .705, -.314,  .155,
#'                    .820, -.377, -.104,
#'                    .661,  .397,  .077,
#'                    .457,  .294, -.488,
#'                    .765,  .428,  .009),
#'                  nrow = 9, ncol = 3, byrow = TRUE)
#'
#' ## Create partially-specified target matrix
#' Targ <- matrix(c(NA, 0,  NA,
#'                  NA, 0,  0,
#'                  NA, 0,  0,
#'                  NA, NA, NA,
#'                  NA, NA, 0,
#'                  NA, NA, NA,
#'                  .7, NA, NA,
#'                  0,  NA, NA,
#'                  .7, NA, NA),
#'                nrow = 9, ncol = 3, byrow = TRUE)
#'
#' ## Perform target rotation
#'                rotate       = "targetT",
#'                targetMatrix = Targ,
#'                digits       = 3)$loadings #' Out3 #' @import GPArotation #' @import stats #' @export faMain <- function(X = NULL, ## Raw data matrix or data.frame R = NULL, ## Correlation matrix for extraction n = NULL, ## Sample Size for faml numFactors = NULL, ## Number of factors to extract facMethod = "fals", ## Factor extraction method urLoadings = NULL, ## Unrotated loading matrix rotate = "oblimin", ## Which rotation to use targetMatrix = NULL, ## If targetT or targetQ, matrix of specified loadings bootstrapSE = FALSE, ## Whether to bootstrap stand. errors numBoot = 1000, ## Number of bootstrapped samples drawn CILevel = .95, ## Bootstrap SE confidence level Seed = 1, ## Seed for bootstrap reproducibility digits = NULL, ## Round all output faControl = NULL, ## Control factor extraction parameters rotateControl = NULL, ## Control rotation tuning parameters ...) { ## ~~~~~~~~~~~~~~~~~~ ## #### Error checking #### ## ~~~~~~~~~~~~~~~~~~ ## if( !is.null(R) ) { # Check if only row or col names given # then remove names rowOrCol <- which.max(c(length(rownames(R)), length(colnames(R)))) Vnames <- dimnames(R)[[rowOrCol]] #Put row and col names back on matrix dimnames(R) <- list(Vnames,Vnames) }# END if(!is.null(R)) FLAGurLoadings <- FALSE if(!is.null(urLoadings)) FLAGurLoadings <- TRUE ## Check rotate ## Is a correct rotation is specified if ( !is.character(rotate) ) { stop("The 'rotate' argument must be a character string. See ?faMain for available options.") } # END if ( !is.character(rotate) ) ## Character string of all available rotate options PlausibleRotations <- c("oblimin", "quartimin", "targetT", "targetQ", "oblimax", "entropy", "quartimax", "varimax", "simplimax", "bentlerT", "bentlerQ", "tandemI", "tandemII", "geominT", "geominQ", "promaxQ", "cfT", "cfQ", "infomaxT", "infomaxQ", "mccammon", "bifactorT", "bifactorQ", "none") ## If rotate argument is incorrectly specified, return error if ( rotate %in% PlausibleRotations == FALSE ) { stop("The 'rotate' argument is incorrectly specified. Check for spelling errors (case sensitive).") } # END if (rotate %in% PlausibleRotations == FALSE) { ## Create nVar object, used in defining 'k' param in cnRotate list if (!is.null(X)) nVar <- ncol(X) if (!is.null(R)) nVar <- ncol(R) if (!is.null(urLoadings)) nVar <- nrow(urLoadings) ## Check R ## Make sure either corr mat or factor structure is supplied if ( is.null(R) & is.null(urLoadings) & is.null(X) ) { stop("Must specify 'X', 'R', or 'urLoadings' arguments.") } # END if ( is.null(R) & is.null(urLoadings) & is.null(X) ) ## If factor analysis is run, must specify numFactors if ( is.null(urLoadings) & is.null(numFactors) ) { stop("The user needs to specify the number of factors to extract.") } # END if ( is.null(urLoadings) & is.null(numFactors) ) ## Check X (if specified) ## Only check the following if X is specified if ( !is.null(X) ) { ## Missingness in data? if (nrow(X) != nrow(X[complete.cases(X),])) { warning("There are missing values in the data frame. See details for how missing data is handled in the cor function.") } # END if (nrow(X) != nrow(X[complete.cases(X),])) ## Class must be matrix or DF to analyze via cor() function if ( !any( class(X) %in% c("matrix", "data.frame", "loadings") ) ){ stop("'X' must be of class matrix, data.frame, or loadings.") } # END if ( class(X) %in% c("matrix", "data.frame", "loadings")) } # END if ( !is.null(X) ) ## Digits ## If a digits argument is not supplied, use R's default option if ( is.null(digits) ) { digits <- options()$digits
} # END if ( is.null(digits) )

## BootstrapSE

## Must give the raw data to perform bootstraps
if ( bootstrapSE == TRUE && is.null(X)) {
stop("The raw data are required to compute the bootstrapped standard errors.")
} # END if ( bootstrapSE == TRUE && is.null(X))

## rotateControl

#### ------- DEFINE cnRotate -------- ####

## Assign the default values for the rotateControl list
cnRotate <- list(numberStarts = 10,
gamma        = 0,
delta        = .01,
kappa        = 0,
k            = nVar,
standardize  = "none",
epsilon      = 1e-5,
power        = 4,
maxItr       = 15000)

## Test that rotateControl arguments are correctly specified
if ( !is.null(rotateControl) ) {

## If rotateControl is specified & standardize is misspecified:
if ( !is.null(rotateControl$standardize) && rotateControl$standardize %in% c("none", "Kaiser", "CM") == FALSE) {

## Produce an error
stop("The 'standardize' argument must be: 'none', 'Kaiser', or 'CM'.")
} # END

## Number of names in the default rotateControl list
cnLength <- length( names(cnRotate) )

## Number of all unique names across user rotateControl and default rotateControl lists
allLength <- length( unique( c( names(rotateControl), names(cnRotate) ) ) )

## If the lengths differ, it is because of spelling issue
if (cnLength != allLength) {

## Find the incorrect arg/s
incorrectArgs <- which( (names(rotateControl) %in% names(cnRotate) ) == FALSE)

# stop("The following arguments are not valid inputs for the list of rotateControl arguments: ", paste0( paste(c(names(rotateControl)[incorrectArgs]), collapse = ", "), "." ) )
stop(paste("The following arguments are not valid inputs for the list of rotateControl arguments:", paste(c(names(rotateControl)[incorrectArgs]), collapse = ", "), collapse = ": " ) )

} # END if (cnLength != allLength)

} # END if ( is.null(rotateControl) )

## Change the default values based on user-specified rotateControl arguments
cnRotate[names(rotateControl)] <- rotateControl

## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ##
#### Define Utility Function(s) ####
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ##

## Generate a random orthonormal starting matrix
randStart <- function(dimension) {
qr.Q(qr(matrix(rnorm(dimension^2), dimension, dimension )))
} # END randStart <- function(dimension) {

## Do the rotations, will be called multiple times
internalRotate <- function(lambda,
rotation,
spinMatrix,
rotateControl) {
## Purpose: Simple rotation wrapper
##
##       rotation:      (character) which rotation to use
##       rotateControl: (list) tuning parameters to pass to rotation
##

## Determine column dimensions
matrixDim <- ncol(lambda)

# ## Perform rotation with the specified parameters
# ## if numberStarts = 1 then rotate from the unrotated
# ## solution
# if (rotateControl$numberStarts == 1) { # RandomSpinMatrix <- diag(matrixDim) # } # END if (rotateControl$numberStarts == 1)
#
# if (rotateControl$numberStarts > 1) { # RandomSpinMatrix <- randStart(matrixDim) # } # END if (rotateControl$numberStarts > 1)

switch(rotation,
"none" = {
Phi      = diag( matrixDim ) ) ## Identity matrix
},
"oblimin" = {
GPArotation::oblimin(lambda,
Tmat      = spinMatrix,
gam       = cnRotate$gamma, normalize = FALSE, eps = cnRotate$epsilon,
maxit     = cnRotate$maxItr) }, "quartimin" = { GPArotation::quartimin(lambda, Tmat = spinMatrix, maxit = cnRotate$maxItr,
eps       = cnRotate$epsilon, normalize = FALSE) }, "targetT" = { GPArotation::targetT(lambda, Tmat = spinMatrix, Target = targetMatrix, normalize = FALSE, eps = cnRotate$epsilon,
maxit     = cnRotate$maxItr) }, "targetQ" = { GPArotation::targetQ(lambda, Tmat = spinMatrix, Target = targetMatrix, normalize = FALSE, eps = cnRotate$epsilon,
maxit     = cnRotate$maxItr) }, "oblimax" = { GPArotation::oblimax(lambda, Tmat = spinMatrix, normalize = FALSE, eps = cnRotate$epsilon,
maxit     = cnRotate$maxItr) }, "entropy" = { GPArotation::entropy(lambda, Tmat = spinMatrix, normalize = FALSE, eps = cnRotate$epsilon,
maxit     = cnRotate$maxItr) }, "quartimax" = { GPArotation::quartimax(lambda, Tmat = spinMatrix, normalize = FALSE, eps = cnRotate$epsilon,
maxit     = cnRotate$maxItr) }, "varimax" = { GPArotation::Varimax(lambda, Tmat = spinMatrix, normalize = FALSE, eps = cnRotate$epsilon,
maxit     = cnRotate$maxItr) }, "simplimax" = { GPArotation::simplimax(lambda, Tmat = spinMatrix, normalize = FALSE, eps = cnRotate$epsilon,
k         = cnRotate$k, maxit = cnRotate$maxItr)
},
"bentlerT" = {
GPArotation::bentlerT(lambda,
Tmat      = spinMatrix,
maxit     = cnRotate$maxItr, eps = cnRotate$epsilon,
normalize = FALSE)
},
"bentlerQ" = {
GPArotation::bentlerQ(lambda,
Tmat      = spinMatrix,
maxit     = cnRotate$maxItr, eps = cnRotate$epsilon,
normalize = FALSE)
},
"tandemI" = {
GPArotation::tandemI(lambda,
Tmat      = spinMatrix,
maxit     = cnRotate$maxItr, eps = cnRotate$epsilon,
normalize = FALSE)
},
"tandemII" = {
GPArotation::tandemII(lambda,
Tmat      = spinMatrix,
maxit     = cnRotate$maxItr, eps = cnRotate$epsilon,
normalize = FALSE)
},
"geominT" = {
GPArotation::geominT(lambda,
Tmat      = spinMatrix,
delta     = cnRotate$delta, normalize = FALSE, eps = cnRotate$epsilon,
maxit     = cnRotate$maxItr) }, "geominQ" = { GPArotation::geominQ(lambda, Tmat = spinMatrix, delta = cnRotate$delta,
normalize = FALSE,
eps       = cnRotate$epsilon, maxit = cnRotate$maxItr)
},
"promaxQ" = {
power       = cnRotate$power, standardize = cnRotate$standardize,
epsilon     = cnRotate$epsilon, maxItr = cnRotate$maxItr)
},
"cfT" = {
GPArotation::cfT(lambda,
Tmat      = spinMatrix,
kappa     = cnRotate$kappa, maxit = cnRotate$maxItr,
eps       = cnRotate$epsilon, normalize = FALSE) }, "cfQ" = { GPArotation::cfQ(lambda, Tmat = spinMatrix, kappa = cnRotate$kappa,
eps       = cnRotate$epsilon, normalize = FALSE, maxit = cnRotate$maxItr)
},
"infomaxT" = {
GPArotation::infomaxT(lambda,
Tmat      = spinMatrix,
normalize = FALSE,
eps       = cnRotate$epsilon, maxit = cnRotate$maxItr)
},
"infomaxQ" = {
GPArotation::infomaxQ(lambda,
Tmat      = spinMatrix,
normalize = FALSE,
eps       = cnRotate$epsilon, maxit = cnRotate$maxItr)
},
"mccammon" = {
GPArotation::mccammon(lambda,
Tmat      = spinMatrix,
normalize = FALSE,
eps       = cnRotate$epsilon, maxit = cnRotate$maxItr)
},
"bifactorT" = {
GPArotation::bifactorT(lambda,
Tmat      = spinMatrix,
normalize = FALSE,
eps       = cnRotate$epsilon, maxit = cnRotate$maxItr)
},
"bifactorQ" = {
GPArotation::bifactorQ(lambda,
Tmat      = spinMatrix,
normalize = FALSE,
eps       = cnRotate$epsilon, maxit = cnRotate$maxItr)
})
} # END internalRotate

## Compute Guttman's factor determinacy indices
GuttmanIndices <- function(Lambda,
PhiMat,
SampCorr) {
## Purpose: Compute Guttman (1955) factor indeterminacy indices
##
##          PhiMat: (Matrix) Factor correlation matrix

## Fator structure (works for either oblique or orthogonal model)
facStruct <- Lambda %*% PhiMat

## Try to invert, if not invertible, gives class 'try-error
Rinv  <-try(solve(SampCorr), silent = TRUE)
## If non-invertible, return NAs instead of returning an error
if ( any( class(Rinv) %in% "try-error") ) {
# warning("\n\nEncountered a singular R matrix when computing factor indeterminancy values\n")
return( rep(NA, ncol(Lambda) ) )
} # END if ( any( class(Rinv) %in% "try-error") )
## Factor indeterminacy solution
FI <- sqrt( diag( t(facStruct) %*% Rinv %*% facStruct))
if(max(FI)>1) FI <- rep(NA, length(FI))
FI
} # END GuttmanIndices

#### ------- VARIABLE NAME RETENTION -------- ####

## Extract names of the indicators to rename final output
varNames <- NULL

## If variable names are provided, retain them
if ( !is.null(X) )          varNames <- colnames(X)
if ( !is.null(R) )          varNames <- colnames(R)

## If variable names is NULL, assign them as var 1 to var nItems
if ( is.null(varNames) ) {

## Specify the number of indicators based on the function inputs
if ( !is.null(X) )         nItem <- ncol(X)
if ( !is.null(R) )         nItem <- ncol(R)

varNames <- paste0("var", seq_len(nItem))
} # END if ( is.null(varNames) )

#### ------- BEGIN FUNCTION -------- ####

## Return a copy of the call function.
CALL <- match.call()

## Compute the correlation matrix
if ( !is.null(X) ) {

## Create correlation matrix (default)
R <- cor(X, ...)

} # END if ( is.null(R) )

#### ------- FACTOR EXTRACTION -------- ####

## September 3, 2019  Check Changes
# needed if no factor extraction occurs
faOut <- vector(mode="list")
faOut$faControl <- NULL ## No model fit stuff faModelFit <- NULL ## Compute communalities before rotation faXh2 <- apply(urLoadings^2, 1, sum) } # END if ( !is.null(urLoadings) ) ## If urLoadings is not supplied, extract it if ( is.null(urLoadings) ) { ## Call faX function to extract factor structure matrix ## September 3, 2019 Check Changes faOut <- vector(mode="list") faOut$faControl <- NULL

faOut <- faX(R          = R,
n          = n,
numFactors = numFactors,
facMethod  = facMethod,
faControl  = faControl)

urLoadings <- faOut$loadings[] ## Communalities from factor extraction faXh2 <- faOut$h2

## Extract all other output

convergedPos <- which(names(faOut$faFit) %in% "converged") names(faOut$faFit)[convergedPos] <- "convergedX"
faModelFit <- faOut$faFit } # END if ( is.null(urLoadings) ) ## If urLoadings is provided and numFactors isn't, specify numFactors if ( is.null(numFactors) && !is.null(urLoadings) ) { numFactors <- ncol(urLoadings) } # END if ( is.null(numFactors) && !is.null(urLoadings) ) #### ------- STANDARDIZATION -------- #### ## Call standardize function Stnd <- faStandardize(method = cnRotate$standardize,

## Extract DvInv for later unstandardization
lambda <- Stnd$lambda #### ------- ROTATION -------- #### ## If only 1 factor, rotate must equal "none" if ( numFactors == 1 ) rotate <- "none" ## Pre-allocate a list for the different attempts starts <- vector("list", cnRotate$numberStarts)

## Set the seed for reproducibility in random spin matrices
set.seed(seed = Seed)

## Create a list of matrices to randomly spin the factor structure
starts <- lapply(starts, function(x) randStart(dimension = ncol(lambda)))

## First start is always from the unrotated factor orientation
starts[[1]] <- diag(ncol(lambda))

## For each list element, do the specified rotate and save all output
starts <- lapply(starts, function(x) {

## Rotate the standardized factor structure matrix
rotatedLambda <- internalRotate(lambda        = lambda,
rotation      = rotate,
spinMatrix    = x,
rotateControl = cnRotate)

## If an orthogonal model, make Phi identity matrix
if ( is.null(rotatedLambda$Phi) ) { ## Identity matrix for Phi rotatedLambda$Phi <- diag(numFactors)

} # END if ( is.null(rotatedLambda$Phi) ) ## Unstandardize the estimated loadings rotatedLambda$loadings[] <- Stnd$DvInv %*% rotatedLambda$loadings[]

rotatedLambda

}) # END lapply(starts, function(x) )

#### ------- ORDER UNIQUE SOLUTIONS -------- ####

## Save the minimized Complexity function for each attempt
##    $Table[,2] is the value of criterion at each iteration, which.min ## finds the minimum criterion value across each attempt ## Evaluate the complexity functions to find smallest value if (rotate != "none") { ## For each random start, find the evaluated Complexity function ComplexityFunc <- sapply(starts, function(attempt) min(attempt$Table[, 2]))

} # if (rotate != "none")

## No Complexity function to evaluate when rotate = "none"
if (rotate == "none") {
ComplexityFunc <- rep(1, cnRotate$numberStarts) } # END if (rotate == "none") { ## Find the sort order (minimum first) of the rotation ## complexity fnc values ## Order func creates an ordered vector with scalars representing ## pre-ordered location ## E.g., if the smallest value is originally the 3rd element, ## order will create a vector with the 3rd element as "1" sortedComplexityOrder <- order(ComplexityFunc, decreasing = FALSE) ## UnSpun is *always* 1st rotation. Find where "1" is in ordered vector UnSpunPosition <- which.min(sortedComplexityOrder) ## Create new list to hold the rotated factor output ## (loadings, phi, complexity fnc etc) uniqueSolutions <- vector("list", cnRotate$numberStarts)

## Create list of output from local solutions
for ( iternumStart in 1:cnRotate$numberStarts ){ ## Determine which element of sortedComplexityOrder to grab ## Start with lowest Complexity value, end with highest num <- sortedComplexityOrder[iternumStart] ## Extract the relevant factor loadings and Phi matrices ## "starts" is the unsorted list of rotated output selected <- starts[[num]] #### ------- FACTOR SORT -------- #### ## June 27, 2019: Do not order factors when Procrustes rotation chosen ## ## sort loadings and Phi matrix for non Procrstes if(rotate !="none"){ if( rotate != "targetT" & rotate != "targetQ"){ sortedSols <- orderFactors(Lambda = selected$loadings,
PhiMat  = selected$Phi, salient = .01, reflect = TRUE) ## Overwrite "selected" list selected$loadings <- sortedSols$Lambda selected$Phi      <- sortedSols$PhiMat }#END if( rotate != "targetT" }#END if(rotate !="none") ## Select the factor loadings uniqueSolutions[[iternumStart]]$loadings <- selected$loadings ## Select factor correlations uniqueSolutions[[iternumStart]]$Phi <- selected$Phi ## Select complexity values uniqueSolutions[[iternumStart]]$RotationComplexityValue <- ComplexityFunc[num]

#### ----- FACTOR INDETERMINACY ----- ####
## compute factor indeterminancy values

## Guttman's factor indeterminacy indices
uniqueSolutions[[iternumStart]]$facIndeterminacy <- GuttmanIndices(Lambda = selected$loadings,
PhiMat = selected$Phi, SampCorr = R ) }else{ uniqueSolutions[[iternumStart]]$facIndeterminacy <- rep(NA, numFactors)
}

# END if(!is.null((R)))

## Did the local optima solution converge
uniqueSolutions[[iternumStart]]$RotationConverged <- selected$convergence

} # END for (iternumStart in 1:length(sortedDiscOrder))

## Extracted rotation complexity values from uniqueSolutions to find solution sets
DisVal <- unlist(lapply(uniqueSolutions, function(x) x$RotationComplexityValue)) ## Round the discrepancy values to specified number of digits DisVal <- round(x = DisVal, digits = nchar(cnRotate$epsilon))

## Determine the sets of local solutions
localMins <- unique(DisVal)

# Classify solutions into UNIQUE groups
nGrp <- length(localMins)
localGrps <- vector("list", nGrp)
names(localGrps) <- as.character(localMins)

## Sort the discrep functions into the local minima buckets
for (local in 1:nGrp) {

## Determine which local solution belong in this set
whichBelong <- which(DisVal == localMins[local])

## Select those matrices, put them in localGrps
localGrps[[local]] <- uniqueSolutions[whichBelong]

} # END for (local in 1:nGrp)

## uniqueSolutions is ordered so 1 is minimum discrep value
minLambda  <- uniqueSolutions[[1]]$loadings minPhi <- uniqueSolutions[[1]]$Phi
## CG EDITS (30 sept 19): Changed "facIndeter" to a column vector
facIndeter <- data.frame("FI" = uniqueSolutions[[1]]facIndeterminacy) ## If factor indeter. not computed, give NA values if ( is.null(facIndeter) ) facIndeter <- rep(NA, ncol(minLambda)) ## If target matrix given then align final factor solution ## to the target if ( !is.null(targetMatrix) && rotate %in% c("targetT", "targetQ", "none") == FALSE ) { minAlign <- faAlign(F1 = targetMatrix, F2 = minLambda, Phi2 = minPhi, MatchMethod = "LS") minLambda <- minAlignF2
minPhi     <- minAlign$Phi2 facIndeter <- data.frame("FI" = facIndeter[ minAlign$FactorMap[2, ], ])
}

## CG EDITS (30 sept 19): Moved communality computation to before bootstraps
## Compute communalities. If orthogonal model, Phi is identity matrix
communalityDF <-
data.frame("h2" = diag(minLambda %*% minPhi %*% t(minLambda)))

#### -------- BOOTSTRAP SETUP -------- ####

## If true, compute bootstrap standard errors
if (bootstrapSE == TRUE) {

# Number of subjects
nSubj <- nrow(X)

# Lists/variables to hold output
fList <-
## Hold factor correlations from bootstraps
phiList <-
## Hold factor indeterminacy indices from bootstraps
FIList <-
## Hold communality estimates from bootstraps
h2List <- vector("list", numBoot)

## Number of rows, used below for bootstrap sampling
rows <- 1:nSubj

#### ----____BOOTSTRAP FOR LOOP -------- ####

## Analyses on 'numBoot' number of random resamples of X matrix
for (iSample in seq_len(numBoot)) {

## Set the seed for reproducibility
set.seed(iSample + Seed)

## Resample (with replacement) from X raw data matrix
bsSample <- sample(x       = rows,
size    = nSubj,
replace = TRUE)

## Create correlation matrix from the bootstrap sample
# Rsamp <- cor(X[bsSample, ], ...)
# if missing data present in the original sample, use
# missing data method passed to faMain fnc

Rsamp <- cor(X[bsSample, ], ...)

## Extract unrotated factors using resampled data matrix
bsLambda <- faX(R          = Rsamp,
numFactors = numFactors,
facMethod  = facMethod,
faControl  = faControl)$loadings ## Conduct standardization bsStnd <- faStandardize(method = cnRotate$standardize,
lambda = bsLambda)

bsLambda <- bsStnd$lambda ## Find the "global" min out of all the random starting configurations bsStarts <- vector("list", cnRotate$numberStarts)

bsStarts <- lapply(bsStarts, function(x) randStart(dimension = numFactors))

## Conduct rotations from random start values
bsStarts <- lapply(bsStarts, function(x) {

## Rotate the bootstrapped samples
bsRotated <- internalRotate(lambda        = bsLambda,
rotation      = rotate,
spinMatrix    = x,
rotateControl = cnRotate)

## If an orthogonal model, turn Phi into an identity matrix
if ( is.null(bsRotated$Phi) ) { ## Identity matrix bsRotated$Phi <- diag(numFactors)

} # END if ( is.null(bsRotated$Phi) ) ## Return all output, not just Phi matrix bsRotated }) # END bsStarts <- lapply(bsStarts, function(x) ## Find minimum of bsSolutions ## Evaluate the minimum disc functions to find smallest value if ( rotate != "none" ) { ## For each random start, find the evaluated discrepancy function bootstrapComplexityFunc <- sapply(bsStarts, function(attempt) min(attempt$Table[, 2]))

}# END if ( rotate != "none" )

## If no rotation, no discrepancy value to evaluate
if (rotate == "none") bootstrapComplexityFunc <- rep(1, cnRotate$numberStarts) ## Of all random configs, determinine which has the lowest criterion value bsLambda <- bsStarts[[which.min(bootstrapComplexityFunc)]]$loadings
bsPhi    <- bsStarts[[which.min(bootstrapComplexityFunc)]]$Phi ## Unstandardize the minimum rotated solution for this bootstrap sample bsLambda <- bsStnd$DvInv %*% bsLambda

## Pre-allocate list
Aligned <- list()

## Factor alignment does not work on 1 factor (and doesn't make sense)
if (numFactors > 1) {
## Align bootstrap sample with global minimum found earlier
Aligned <- faAlign(F1   = minLambda,
F2   = bsLambda,
Phi2 = bsPhi)

## Extract aligned elements
AlignedLambda <- Aligned$F2 AlignedPhi <- Aligned$Phi2
} # END if (numFactors > 1)

## Cannot align 1-factor model, but can reflect the factor
if (numFactors == 1) {
if ( sum(bsLambda) < 0 ) bsLambda <- bsLambda * -1

## Define newly aligned elements
AlignedLambda <- bsLambda
AlignedPhi    <- bsPhi
} # END if (numFactors == 1)

fList[[iSample]]   <- AlignedLambda

## Save correlations as the bootstrapped factor correlations
phiList[[iSample]] <- AlignedPhi

## Save factor indeterminacy indices
FIList[[iSample]] <- GuttmanIndices(Lambda = AlignedLambda,
PhiMat = AlignedPhi,
SampCorr = Rsamp)

## CG EDITS (30 sept 19): Added list of computed h^2 values
## Save communality estimates
h2List[[iSample]] <- diag(AlignedLambda %*% AlignedPhi %*% t(AlignedLambda))

} # END for (iSample in seq_len(numBoot))

#### -----____ BootStrap STANDARD ERRORS ----- ####

# Convert list of matrices into array
loadArray <- array(unlist(fList), c(nVar, numFactors, numBoot))
phiArray  <- array(unlist(phiList), c(numFactors, numFactors, numBoot))
FIArray   <- array(unlist(FIList), c(1, numFactors, numBoot))
## CG EDITS (30 sept 19): Added array for h2 values
h2Array   <- array(unlist(h2List), c(nVar, 1, numBoot))

# Bootstrap standard errors for factor correlations
phiSE <- apply(phiArray, 1:2, sd)
phiSE <- round(phiSE, digits)

## CG EDITS (30 sept 19): Added bootstrap SEs for communalities
h2SE <- apply(h2Array, 1, sd)

## ----- CONFIDENCE INTERVALS ----- ##

# Set alpha level
alpha <- 1 - CILevel

# Confidence interval matrices for factor correlations
phiCI.upper <- apply(phiArray, 1:2, quantile, 1 - (alpha / 2))
phiCI.lower <- apply(phiArray, 1:2, quantile, (alpha / 2))

## Round Phi
phiCI.upper <- round(phiCI.upper, digits)
phiCI.lower <- round(phiCI.lower, digits)

## CG EDITS (30 sept 19): Added CIs for fact indeterminacy estimates
## ----____Bootstrap SE FI----
if(!any(is.na(FIArray))){
FISE <- apply(FIArray, 1:2, sd)
# FISE <- round(FISE, digits)
FICI.upper <- apply(FIArray, 2, quantile, 1 - (alpha / 2))
FICI.lower <- apply(FIArray, 2, quantile, (alpha / 2))

## CG EDITS (30 sept 19): Create data frame of fac indeterminacy estimates
facIndeter <- data.frame(facIndeter,
t(FISE),
FICI.lower,
FICI.upper)
## CG EDITS (30 sept 19): Add column names to data frame of FI estimates
colnames(facIndeter) <- c("FI",
"SE",
paste0((alpha / 2) * 100, "th percentile"),
paste0((1 - (alpha / 2)) * 100, "th percentile"))
rownames(facIndeter) <-  paste0("f", 1:numFactors)
}else{
facIndeter <- NA
FISE <- NA
}

## CG EDITS (30 sept 19): Added CIs for communality estimates
h2CI.upper <- apply(h2Array, 1, quantile, 1 - (alpha / 2))
h2CI.lower <- apply(h2Array, 1, quantile, (alpha / 2))

## CG EDITS (30 sept 19): Create data frame of communalities
communalityDF <- data.frame(communalityDF, ## Obtained communalities
h2SE,          ## Communality stand. errors
h2CI.lower,    ## Lower CI bound
h2CI.upper)    ## Upper CI bound

## CG EDITS (30 sept 19): Add column names to data frame of communalities
colnames(communalityDF) <-
c("h2",
"SE",
paste0((alpha / 2) * 100, "th percentile"),
paste0((1 - (alpha / 2)) * 100, "th percentile"))

#### ------- NAME BOOTSTRAP OBJECTS -------- ####

## Dimension names

## Name the factor indicators

## Name the factors
colnames(phiSE) <- rownames(phiSE) <-
colnames(phiCI.upper) <- rownames(phiCI.upper) <-
colnames(phiCI.lower) <- rownames(phiCI.lower) <-
paste0("f", 1:numFactors)

} # END if (bootstrapSE == TRUE)

## If not specified, return NULL
if (bootstrapSE == FALSE) {
phiSE         <-
phiCI.upper   <-
phiCI.lower   <-
phiArray      <-
FIArray       <-
FISE          <-
FICI.upper    <-
FICI.lower    <-
h2SE          <-
h2CI.upper    <-
h2CI.lower    <- NULL
} # END if (bootstrapSE == FALSE)

#### ------- BOOTSTRAP COMMUNALITIES -------- ####

## Check to ensure communalityDF and faH2 are equivalent
CheckEquiv <- all.equal(communalityDF$h2, faXh2, check.attributes = FALSE, ## Ignoring attribute names tolerance = 1e-6) ## precision for equivalence ## If unequivalent, give an ominous warning. if (CheckEquiv == FALSE) { warning("Results may be untrustworthy. Unrotated communalities differ from rotated communalities.") } # END if (CheckEquiv == FALSE) #### ------- NAME REMAINING OBJECTS -------- #### ## Designate the factor names facNames <- paste0("f", 1:numFactors) ## Name objects to have factor indicator names rownames(minLambda) <- ## CG EDITS (30 sept 19): Add row names to data frame of communalities rownames(communalityDF) <- varNames ## Name objects to have factor names ## Change made March 4, 2019 NGW ## handle 1-fac models if(numFactors > 1) colnames(minLambda) <- facNames if(numFactors == 1) dimnames(minLambda) <- list(varNames, facNames) rownames(minPhi) <- colnames(minPhi) <- facNames ## facIndeter is NA if urLoadings is specified if ( any(is.na(facIndeter)) == FALSE) rownames(facIndeter) <- facNames ## Add rotation convergence status to faFit faModelFit$convergedR <- uniqueSolutions[[1]]$RotationConverged #### ----- OUTPUT ----- #### fout <- list(R = R, loadings = minLambda, Phi = minPhi, facIndeterminacy = facIndeter, h2 = communalityDF, loadingsSE = loadSE, CILevel = CILevel, loadingsCIupper = loadCI.upper, loadingsCIlower = loadCI.lower, PhiSE = phiSE, PhiCIupper = phiCI.upper, PhiCIlower = phiCI.lower, facIndeterminacySE = FISE, localSolutions = uniqueSolutions, numLocalSets = nGrp, localSolutionSets = localGrps, loadingsArray = loadArray, PhiArray = phiArray, facIndeterminacyArray = FIArray, faControl = faOut$faControl,
faFit                 = faModelFit,
rotate                = rotate,
rotateControl         = cnRotate,
unSpunSolution        = uniqueSolutions[[UnSpunPosition]],
targetMatrix          = targetMatrix,
Call                  = CALL)

class(fout) <- 'faMain'
fout

} ## END faMAIN


## Try the fungible package in your browser

Any scripts or data that you put into this service are public.

fungible documentation built on Sept. 29, 2021, 1:06 a.m.