Nothing
R/faIB.R
In fungible: Psychometric Functions from the Waller Lab
Defines functions
faIB
Documented in
faIB
#' Inter-Battery Factor Analysis by the Method of Maximum Likelihood
#'
#' This function conducts maximum likelihood inter-battery factor analysis using procedures described by Browne (1979).
#' The unrotated solution can be rotated (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) structured in a subject
#' (row) by variable (column) format. Defaults to \code{X = NULL}.
#' @param R (Matrix) A correlation matrix. Defaults to \code{R = NULL}.
#' @param n (Numeric) Sample size associated with either the raw data (X) or
#' the correlation matrix (R). Defaults to \code{n = NULL}.
#' @param NVarX (Integer) Given batteries X and Y, \code{NVarX} denotes the number of variables in battery X.
#' @param numFactors (Numeric) The number of factors to extract for subsequent
#' rotation. Defaults to \code{numFactors = NULL}.
#' @param itemSort (Logical) if \code{itemSort = TRUE} the factor loadings will be sorted within batteries.
#' @param rotate (Character) Designate which rotation algorithm to apply. The
#' following are available rotation options: "oblimin", "quartimin",
#' "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 Seed (Integer) Starting seed for the random number generator.
#' @inheritParams faMain
#'
#' @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{faIB} function will produce abundant output in addition
#' to the rotated inter-battery factor pattern and factor correlation matrices.
#' \itemize{
#' \item \strong{loadings}: (Matrix) The rotated inter-battery factor solution with the
#' 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{fit}: (Vector) A vector containing the following fit statistics:
#' \itemize{
#' \item \strong{chiSq}: Chi-square goodness of fit value (see Browne, 1979, for details). Note that we apply Lawley's (1959) correction when computing the chi-square value.
#' \item \strong{DF}: Degrees of freedom for the estimated model.
#' \item \strong{p-value}: P-value associated with the above chi-square statistic.
#' \item \strong{MAD}: Mean absolute difference between the model-implied and the sample across-battery correlation matrices. A lower value indicates better fit.
#' \item \strong{AIC}: Akaike's Information Criterion where a lower value indicates better fit.
#' \item \strong{BIC}: Bayesian Information Criterion where a lower value indicates better fit.
#' }
#' \item \strong{R}: (Matrix) Returns the (possibly sorted) correlation matrix, useful when raw data are supplied.
#' If \code{itemSort = TRUE} then the returned matrix is sorted to be consistent with the factor loading matrix.
#' \item \strong{Rhat}: (Matrix) The (possibly sorted) reproduced correlation matrix.If \code{itemSort = TRUE} then the returned matrix is sorted to be consistent with the factor loading matrix.
#' \item \strong{Resid}: (Matrix) A (possibly sorted) residual matrix (R - Rhat) for the between battery correlations.
#' \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{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{loadings}: (Matrix) the factor loadings,
#' \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.
#' \item \strong{rotate} (Character) The chosen rotation algorithm.
#' \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{Call}: (call) A copy of the function call.
#' }
#'
#' @references Boruch, R. F., Larkin, J. D., Wolins, L., & MacKinney, A. C. (1970).
#' Alternative methods of analysis: Multitrait-multimethod data. \emph{Educational
#' and Psychological Measurement, 30}(4), 833–853.
#' https://doi.org/10.1177/0013164470030004055
#' @references Browne, M. W. (1979). The maximum-likelihood solution in inter-battery factor analysis.
#' \emph{British Journal of Mathematical and Statistical Psychology, 32}(1), 75-86.
#' @references Browne, M. W. (1980). Factor analysis of multiple batteries by maximum likelihood.
#' \emph{British Journal of Mathematical and Statistical Psychology, 33}(2), 184-199.
#' @references Browne, M. W. (2001). An overview of analytic rotation in
#' exploratory factor analysis. \emph{Multivariate Behavioral Research, 36}(1), 111-150.
#' @references Burnham, K. P. & Anderson, D. R. (2004). Multimodel inference: Understanding AIC and BIC in model selection.
#' \emph{Sociological methods and research, 33}, 261-304.
#' @references Cudeck, R. (1982). Methods for estimating between-battery factors,
#' \emph{Multivariate Behavioral Research, 17}(1), 47-68. 10.1207/s15327906mbr1701_3
#' @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 Tucker, L. R. (1958). An inter-battery method of factor analysis.
#' \emph{Psychometrika, 23}(2), 111-136.
#' @family Factor Analysis Routines
#'
#' @author
#' \itemize{
#' \item Niels G. Waller (nwaller@umn.edu)
#' \item Casey Giordano (Giord023@umn.edu)
#'}
#'
#' @examples
#'
#' # Example 1:
#' # Example from: Browne, M. W. (1979).
#' #
#' # Data originally reported in:
#' # Thurstone, L. L. & Thurstone, T. G. (1941). Factorial studies
#' # of intelligence. Psychometric Monograph (2), Chicago: Univ.
#' # Chicago Press.
#'
#' R.XY <- matrix(c(
#' 1.00, .554, .227, .189, .461, .506, .408, .280, .241,
#' .554, 1.00, .296, .219, .479, .530, .425, .311, .311,
#' .227, .296, 1.00, .769, .237, .243, .304, .718, .730,
#' .189, .219, .769, 1.00, .212, .226, .291, .681, .661,
#' .461, .479, .237, .212, 1.00, .520, .514, .313, .245,
#' .506, .530, .243, .226, .520, 1.00, .473, .348, .290,
#' .408, .425, .304, .291, .514, .473, 1.00, .374, .306,
#' .280, .311, .718, .681, .313, .348, .374, 1.00, .672,
#' .241, .311, .730, .661, .245, .290, .306, .672, 1.00), 9, 9)
#'
#'
#' dimnames(R.XY) <- list(c( paste0("X", 1:4),
#' paste0("Y", 1:5)),
#' c( paste0("X", 1:4),
#' paste0("Y", 1:5)))
#'
#' out <- faIB(R = R.XY,
#' n = 710,
#' NVarX = 4,
#' numFactors = 2,
#' itemSort = FALSE,
#' rotate = "oblimin",
#' rotateControl = list(standardize = "Kaiser",
#' numberStarts = 10),
#' Seed = 1)
#'
#' # Compare with Browne 1979 Table 2.
#' print(round(out$loadings, 2))
#' cat("\n\n")
#' print(round(out$Phi,2))
#' cat("\n\n MAD = ", round(out$fit["MAD"], 2),"\n\n")
#' print( round(out$facIndeterminacy,2) )
#'
#'
#' # Example 2:
#' ## Correlation values taken from Boruch et al.(1970) Table 2 (p. 838)
#' ## See also, Cudeck (1982) Table 1 (p. 59)
#' corValues <- c(
#' 1.0,
#' .11, 1.0,
#' .61, .47, 1.0,
#' .42, -.02, .18, 1.0,
#' .75, .33, .58, .44, 1.0,
#' .82, .01, .52, .33, .68, 1.0,
#' .77, .32, .64, .37, .80, .65, 1.0,
#' .15, -.02, .04, .08, .12, .11, .13, 1.0,
#' -.04, .22, .26, -.06, .07, -.10, .07, .09, 1.0,
#' .13, .21, .23, .05, .07, .06, .12, .64, .40, 1.0,
#' .01, .04, .01, .16, .05, .07, .05, .41, -.10, .29, 1.0,
#' .27, .13, .18, .17, .27, .27, .27, .68, .18, .47, .33, 1.0,
#' .24, .02, .12, .12, .16, .23, .18, .82, .08, .55, .35, .76, 1.0,
#' .20, .18, .16, .17, .22, .11, .29, .69, .20, .54, .34, .68, .68, 1.0)
#'
#' ## Generate empty correlation matrix
#' BoruchCorr <- matrix(0, nrow = 14, ncol = 14)
#'
#' ## Add upper-triangle correlations
#' BoruchCorr[upper.tri(BoruchCorr, diag = TRUE)] <- corValues
#' BoruchCorr <- BoruchCorr + t(BoruchCorr) - diag(14)
#'
#' ## Add variable names to the correlation matrix
#' varNames <- c("Consideration", "Structure", "Sup.Satisfaction",
#' "Job.Satisfaction", "Gen.Effectiveness", "Hum.Relations", "Leadership")
#'
#' ## Distinguish between rater X and rater Y
#' varNames <- paste0(c(rep("X.", 7), rep("Y.", 7)), varNames)
#'
#' ## Add row/col names to correlation matrix
#' dimnames(BoruchCorr) <- list(varNames, varNames)
#'
#' ## Estimate a model with one, two, and three factors
#' for (jFactors in 1:3) {
#' tempOutput <- faIB(R = BoruchCorr,
#' n = 111,
#' NVarX = 7,
#' numFactors = jFactors,
#' rotate = "oblimin",
#' rotateControl = list(standardize = "Kaiser",
#' numberStarts = 100))
#'
#' cat("\nNumber of inter-battery factors:", jFactors,"\n")
#' print( round(tempOutput$fit,2) )
#' } # END for (jFactors in 1:3)
#'
#' ## Compare output with Cudeck (1982) Table 2 (p. 60)
#' BoruchOutput <-
#' faIB(R = BoruchCorr,
#' n = 111,
#' NVarX = 7,
#' numFactors = 2,
#' rotate = "oblimin",
#' rotateControl = list(standardize = "Kaiser"))
#'
#' ## Print the inter-battery factor loadings
#' print(round(BoruchOutput$loadings, 3))
#' print(round(BoruchOutput$Phi, 3))
#'
#'
#'
#' @import GPArotation
#' @import stats
#' @export
faIB <-
function(X = NULL,
R = NULL,
n = NULL,
NVarX = 4,
numFactors = 2,
itemSort = FALSE,
rotate = "oblimin",
bootstrapSE = FALSE,
numBoot = 1000, ## Number of bootstrapped samples drawn
CILevel = .95, ## Bootstrap SE confidence level
rotateControl = NULL,
Seed = 1){
## Return a copy of the call function.
Call <- match.call()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# ----ERROR CHECKING ----
## 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))
## ____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) {
#----____END ERROR CHECKING----
## Set the seed for reproducibility
set.seed(Seed)
# Calculate R matrix if not given
if(!is.null(X)) R <- cor(X)
nVar <- ncol(R)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ##
#### DEFINE FUNCTIONS ####
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ##
## Function to compute model fit indices
ComputeModelFit <- function(R,
NVarX,
numFactors,
Lambda,
Phi){
## Define the number of variables in the model
nVar <- nrow(Lambda)
## Compute model-implied (reduced) correlation matrix
Rhat <- Lambda %*% Phi %*% t(Lambda)
## Find the residual correlation matrix
Resid <- R - Rhat
## Residual XY correlations
ResidRXY <- Resid[1:NVarX, (NVarX+1):nVar]
## Mean absolute difference of the residual XY correlation matrix
MAD <- mean( abs(ResidRXY) )
## Redefine the number of X variables
p <- NVarX
## Within X-battery correlation matrix
R.XX <- R[1:p, 1:p]
## Within Y-Battery correlation matrix
R.YY <- R[(p+1):nVar, (p+1):nVar]
## Between battery correlation matrix
R.XY <- R[1:p, (p+1):nVar]
# Squared canonical correlations
rhoSq <- Re(eigen(R.XY %*% solve(R.YY) %*% t(R.XY) %*% solve(R.XX))$values)
# Browne 1979 formula for corrected chi-square
chiSqTmp <- 1
for(i in (numFactors+1):NVarX ){
chiSqTmp <- chiSqTmp * (1 - rhoSq[i])
}
# see Browne 1979 p. 84
pstarInv <- 0
for(i in 1:numFactors){
pstarInv <- pstarInv + 1/rhoSq[i]
}
# Lawley 1959 correction for n
## NOTE: No correction if sample size (n) is not specified
if ( !is.null(n) ) {
nstar <- (n - 1 - numFactors - .5*(nVar+1) + pstarInv)
chiSq <- -nstar * log(chiSqTmp)
dfChiSq <- (NVarX - numFactors)*( (nVar - NVarX) - numFactors)
if(dfChiSq <= 0) warning("\n\n*****Zero Degrees of Freedom*****\n")
# check if rotation is orthogonal
# orthogonal is a logical (T/F)
orthogonal <- all.equal(Phi, diag(numFactors),
check.attributes = FALSE,
check.names = FALSE) == TRUE
# ~~~~~~~~~~~~~~AIC and BIC~~~~~~~~~~~~~~~~~
# qk = parsimony corrections factor
if (orthogonal) {
qk <- numEstParams <- (nVar * numFactors)
} else {
qk <- numEstParams <- (nVar * numFactors) + numFactors*(numFactors - 1)/2
} # END if (orthogonal)
#BIC
BIC <- chiSq + log(n) * qk
#AIC
AIC <- chiSq + 2 * qk
} else {
## If sample size is NULL, set the following as NA
BIC <- AIC <- chiSq <- dfChiSq <- NA
} # END if ( !is.null(n) )
## Return
list(Rhat = Rhat,
Resid = Resid,
MAD = MAD,
chiSq = chiSq,
DF = dfChiSq,
AIC = AIC,
BIC = BIC)
}# ED Compute Model Fit
## Extract names of the indicators to rename final output
GetVarNames <- function(){
varNames <- NULL
## If variable names are provided, retain them
if ( !is.null(X) ) varNames <- colnames(X)
if ( !is.null(R) ) varNames <- colnames(R)
if(is.null(varNames )) varNames <- paste0("V", 1:nVar)
varNames
}
# Main code to compute MLE multiple battery factor loadings
MBFact <- function(R, NVarX){
NVar <- ncol(R)
p <- NVarX
R.XX <- R[1:p, 1:p]
R.YY <- R[(p+1):NVar, (p+1):NVar]
R.XY <- R[1:p, (p+1):NVar]
F <- solve(R.XX) %*% R.XY %*% solve(R.YY) %*% t(R.XY)
VDV.F <- eigen(F)
V.F <- VDV.F$vectors
D.F <- VDV.F$values
G <- diag(sqrt(diag(t(V.F) %*% R.XX %*% V.F)))
B.X <- V.F %*% solve(G)
B.Y <- solve(R.YY)%*% t(R.XY) %*% B.X %*% diag(1/sqrt(D.F))
# Unrotated loadings for X (sqrt of sigular values)
LX <-(R.XX %*% B.X)[, 1:numFactors, drop = FALSE] %*%
diag(D.F^.25)[1:numFactors, 1:numFactors]
# Unrotated loadings for Y
LY <- R.YY %*% B.Y[, 1:numFactors, drop = FALSE] %*%
diag(D.F^.25)[1:numFactors, 1:numFactors]
# &&&& Why do we have to ask for the Real part of the complex number
urLoadings <- Re(rbind(LX, LY))
urLoadings
}#END MBFact
SortItems <- function(Lambda1, varNames, NVarX){
soutX <- faSort(Lambda1[1:NVarX, ],
phi = NULL,
BiFactor = FALSE,
salient = 0.25,
reflect = FALSE)
Lambda1[1:NVarX, ] <- soutX$loadings
soutY <- faSort(Lambda1[(NVarX + 1):nVar, ],
phi = NULL,
BiFactor = FALSE,
salient = 0.25,
reflect = FALSE)
Lambda1[(NVarX + 1):nVar, ] <- soutY$loadings
newOrder <- c(soutX$sortOrder,
(soutY$sortOrder + NVarX))
rownames(Lambda1) <- varNames[newOrder]
list( Lambda1 = Lambda1,
newOrder = newOrder)
}# END SortItems
#----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
## Compute Guttman's factor determinacy indices
GuttmanIndices <- function(Lambda,
PhiMat,
SampCorr) {
## Purpose: Compute Guttman (1955) factor indeterminacy indices
##
## Args: Lambda: (Matrix) Rotated factor loadings matrix
## PhiMat: (Matrix) Factor correlation matrix
## Fator structure (works for either oblique or orthogonal model)
facStruct <- Lambda %*% PhiMat
Rinv <-try(solve(SampCorr), silent = TRUE)
## If non-invertible, return NA values 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
## 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
##
## Args: lambda: (matrix) unrotated loadings to rotate
## rotation: (character) which rotation to use
## spinMatrix: (Matrix) randomly spin unrotated factor loadings
## 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" = {
list(loadings = lambda,
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" = {
promaxQ(urLoadings = lambda,
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
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# ---- MAIN PROGRAM----
# ----____Get varNames----
varNames <- GetVarNames()
#----____Compute MB factor loadings----
urLoadings <- MBFact(R, NVarX)
## If only 1 factor, rotate must equal "none"
if ( numFactors == 1 ) rotate <- "none"
## ----____ Standardize loadings ####
Stnd <- faStandardize(method = cnRotate$standardize,
lambda = urLoadings)
## Extract DvInv for later unstandardization
lambda <- Stnd$lambda
##----____ Rotate loadings from random spins ####
## 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(randSpinMatrix) {
rotatedLambda <- internalRotate(lambda = lambda,
rotation = rotate,
spinMatrix = randSpinMatrix,
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[]
## Return whole list, not just 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]]
# December 9, 2019 Note that order factors (Factor Sort) is not yet implimented
#### -------____ Sort rotated factor solutions ####
## 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]
## ----____Compute Factor Indeterminancy ####
## if the user provides urloadings then do not
## compute factor indeterminancy values
## Guttman's factor indeterminacy indices
uniqueSolutions[[iternumStart]]$facIndeterminacy <-
GuttmanIndices(Lambda = selected$loadings,
PhiMat = selected$Phi,
SampCorr = R ) # perhaps change name of R matrix
## 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 create solution sets
DisVal <- unlist(lapply(uniqueSolutions, function(x) x$RotationComplexityValue))
## Determine the number of digits for rotation convergence
## NOTE: subtract 2 because the leading zero and decimal count as characters
## NOTE: Convert scientific notation (e.g., 1e-5) to decimal form (works if non-scientific notation)
numDigits <- nchar( format(cnRotate$epsilon, scientific = FALSE) ) - 2
## Round the discrepancy values to specified number of digits
DisVal <- round(x = DisVal,
digits = numDigits)
## 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)
## Select the loading matrix with the lowest discrep value
## 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("FacIndeterminacy" = uniqueSolutions[[1]]$facIndeterminacy)
## If factor indeter. not computed, give NA values
if ( is.null(facIndeter) ) facIndeter <- rep(NA, ncol(minLambda))
Lambda1 <- minLambda
## ----____Reflect factors ----
# If reflect = TRUE then reflect factors s.t.
# salient loadings are positive
## Determine whether factors are negatively oriented
if (numFactors > 1) {
Dsgn <- diag(sign(colSums(Lambda1^3)))
} else {
Dsgn <- matrix(sign(colSums(Lambda1^3)), 1, 1)
} # END if (numFactors > 1)
## If factors negatively oriented, multiply by -1, else multiply by 1
Lambda1 <- Lambda1 %*% Dsgn
if (!is.null(minPhi)) {
## If factor is negative, reverse corresponding factor correlations
minPhi <- Dsgn %*% minPhi %*% Dsgn
} # END if (!is.null(minPhi))
# ----____Sort Items?----
# if no item sort, newOrder is original order
newOrder <- 1:nVar
if(itemSort == TRUE){
if(numFactors == 1){
newOrderX <- sort.list(Lambda1[1:NVarX])
newOrderY <- sort.list(Lambda1[(NVarX+1):nVar])
newOrder <- c(newOrderX , (newOrderY+NVarX) )
Lambda1 <- matrix(Lambda1[newOrder], nVar, 1)
rownames(Lambda1) <- varNames[newOrder]
colnames(Lambda1) <- "F1"
PhiNames <- colnames(Lambda1) <- paste0("F", 1:ncol(Lambda1))
dimnames(minPhi) <- list(PhiNames, PhiNames)
}
if(numFactors > 1){
sorted_out <- SortItems(Lambda1, varNames, NVarX)
Lambda1 <- sorted_out$Lambda1
newOrder <- sorted_out$newOrder
PhiNames <- colnames(Lambda1) <- paste0("F", 1:ncol(Lambda1))
dimnames(minPhi) <- list(PhiNames, PhiNames)
}
}else{
rownames(Lambda1) <- varNames
PhiNames <- colnames(Lambda1) <- paste0("F", 1:ncol(Lambda1))
dimnames(minPhi) <- list(PhiNames, PhiNames)
} #END Sort Items?
# if itemSort = TRUE the reOrder items in R
R <- R[newOrder, newOrder]
#---- ____Compute Model Fit----
fitOut <- ComputeModelFit(R = R,
NVarX = NVarX,
numFactors = numFactors,
Lambda = Lambda1,
Phi = minPhi)
Rhat <- fitOut$Rhat
Resid <- fitOut$Resid
MAD <- fitOut$MAD
chiSq <- fitOut$chiSq
DF <- fitOut$DF
AIC <- fitOut$AIC
BIC <- fitOut$BIC
#### -------- BOOTSTRAP SETUP -------- ####
## If true, compute bootstrap standard errors
if (bootstrapSE == TRUE) {
#This reorders the cols of X if itemSort = TRUE
X <- X[, newOrder]
# Number of subjects
nSubj <- nrow(X)
# Lists/variables to hold output
## Hold factor loadings from bootstraps
fList <-
## Hold factor correlations from bootstraps
phiList <-
## Hold factor indeterminacy indices from bootstraps
FIList <- 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, ])
## ____ Extract factor loadings ####
## Extract unrotated factors using resampled data matrix
bsLambda <- MBFact(R = Rsamp,
NVarX = NVarX)
## Conduct standardization
bsStnd <- faStandardize(method = cnRotate$standardize,
lambda = bsLambda)
## Extract the standardized bootstrapped (unrotated) factor loadings
bsLambda <- bsStnd$lambda
## ____ Rotate from random spins ####
## Find the "global" min out of all the random starting configurations
## Initialize a list to store rotations from random spins
bsStarts <- vector("list", cnRotate$numberStarts)
## Populate the list with random orthogonal matrices (i.e., start values for rotation)
bsStarts <- lapply(bsStarts, function(x) randStart(dimension = numFactors))
## Conduct rotations from random start values
bsStarts <- lapply(bsStarts, function(randSpinMatrix) {
## Rotate the bootstrapped samples
bsRotated <- internalRotate(lambda = bsLambda,
rotation = rotate,
spinMatrix = randSpinMatrix,
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
bsRotated
}) # END bsStarts <- lapply(bsStarts, function(x)
## ____ Order unique solutions ####
## 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
bsMinimum <- which.min(bootstrapComplexityFunc)
## Extract global minimum factor loadings and Phi matrices from
bsLambda <- bsStarts[[bsMinimum]]$loadings
bsPhi <- bsStarts[[bsMinimum]]$Phi
## Unstandardize the 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 = Lambda1,
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^3) < 0 ) bsLambda <- bsLambda * -1
## Define newly aligned elements
AlignedLambda <- bsLambda
AlignedPhi <- bsPhi
} # END if (numFactors == 1)
## Save loadings as the bootstrapped factor loadings
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)
} # END for (iSample in seq_len(numBoot))
#### -----____ BootStrap stnd 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))
## Bootstrap standard errors for factor loadings
loadSE <- apply(loadArray, 1:2, sd)
# Bootstrap standard errors for factor correlations
phiSE <- apply(phiArray, 1:2, sd)
## ____ Confidence intervals ####
## Set alpha level
alpha <- 1 - CILevel
## Confidence interval matrices for factor loadings
loadCI.upper <- apply(loadArray, 1:2, quantile, 1 - (alpha / 2))
loadCI.lower <- apply(loadArray, 1:2, quantile, (alpha / 2))
# 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))
## ____ Bootstrap Fac Indeterminacy SE ####
## If none of the fac indeter values are NA, proceed. Else, set as NA
if (!any(is.na(FIArray))) {
## Compute the standard error of the bootstrap factor indetermancy values
FISE <- apply(FIArray, 1:2, sd)
## Find the upper/lower bounds of the confidence interval
FICI.upper <- apply(FIArray, 2, quantile, 1 - (alpha / 2))
FICI.lower <- apply(FIArray, 2, quantile, (alpha / 2))
## Create data frame of fac indeterminacy estimates
facIndeter <- data.frame(facIndeter,
t(FISE),
FICI.lower,
FICI.upper)
## 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"))
## Add row names to correspond to each factor
rownames(facIndeter) <- paste0("f", 1:numFactors)
}else{
facIndeter <- NA
FISE <- NA
} # END if (!any(is.na(FIArray)))
#### ------- NAME BOOTSTRAP OBJECTS -------- ####
## Dimension names
## Name the factor indicators
rownames(loadSE) <-
rownames(loadCI.upper) <-
rownames(loadCI.lower) <- varNames[newOrder]
## Name the factors
colnames(loadSE) <-
colnames(loadCI.upper) <-
colnames(loadCI.lower) <-
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) {
loadSE <-
loadCI.upper <-
loadCI.lower <-
loadArray <-
phiSE <-
phiCI.upper <-
phiCI.lower <-
phiArray <-
FIArray <-
FISE <-
FICI.upper <-
FICI.lower <- NULL
} # END if (bootstrapSE == FALSE)
## CG Edits (13 dec 2019): If sample size is NULL, pchisq cannot be computed
if ( !is.null(n) ) {
fitVec <- c(chiSq, DF, pchisq(chiSq, DF), MAD, AIC, BIC)
} else {
fitVec <- c(chiSq, DF, NA, MAD, AIC, BIC)
} # END if ( !is.null(n) )
names(fitVec) <- c("chiSq", "DF", "p-value", "MAD", "AIC", "BIC")
#----____Output ----
invisible(list(loadings = Lambda1, # this is min complexity solutution
Phi = minPhi,
fit = fitVec,
R = R,
Rhat = Rhat,
Resid = Resid,
facIndeterminacy = facIndeter,
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,
rotate = rotate,
loadingsArray = loadArray,
PhiArray = phiArray,
facIndeterminacyArray = FIArray,
rotateControl = cnRotate,
unSpunSolution = uniqueSolutions[[UnSpunPosition]],
Call = Call
))#END inivisble(list( . . .
} ##END faMB
Try the fungible package in your browser
Any scripts or data that you put into this service are public.
fungible documentation built on March 31, 2023, 5:47 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions faIB
Documented in faIB
#' Inter-Battery Factor Analysis by the Method of Maximum Likelihood
#'
#' This function conducts maximum likelihood inter-battery factor analysis using procedures described by Browne (1979).
#' The unrotated solution can be rotated (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) structured in a subject
#' (row) by variable (column) format. Defaults to \code{X = NULL}.
#' @param R (Matrix) A correlation matrix. Defaults to \code{R = NULL}.
#' @param n (Numeric) Sample size associated with either the raw data (X) or
#' the correlation matrix (R). Defaults to \code{n = NULL}.
#' @param NVarX (Integer) Given batteries X and Y, \code{NVarX} denotes the number of variables in battery X.
#' @param numFactors (Numeric) The number of factors to extract for subsequent
#' rotation. Defaults to \code{numFactors = NULL}.
#' @param itemSort (Logical) if \code{itemSort = TRUE} the factor loadings will be sorted within batteries.
#' @param rotate (Character) Designate which rotation algorithm to apply. The
#' following are available rotation options: "oblimin", "quartimin",
#' "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 Seed (Integer) Starting seed for the random number generator.
#' @inheritParams faMain
#'
#' @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{faIB} function will produce abundant output in addition
#' to the rotated inter-battery factor pattern and factor correlation matrices.
#' \itemize{
#' \item \strong{loadings}: (Matrix) The rotated inter-battery factor solution with the
#' 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{fit}: (Vector) A vector containing the following fit statistics:
#' \itemize{
#' \item \strong{chiSq}: Chi-square goodness of fit value (see Browne, 1979, for details). Note that we apply Lawley's (1959) correction when computing the chi-square value.
#' \item \strong{DF}: Degrees of freedom for the estimated model.
#' \item \strong{p-value}: P-value associated with the above chi-square statistic.
#' \item \strong{MAD}: Mean absolute difference between the model-implied and the sample across-battery correlation matrices. A lower value indicates better fit.
#' \item \strong{AIC}: Akaike's Information Criterion where a lower value indicates better fit.
#' \item \strong{BIC}: Bayesian Information Criterion where a lower value indicates better fit.
#' }
#' \item \strong{R}: (Matrix) Returns the (possibly sorted) correlation matrix, useful when raw data are supplied.
#' If \code{itemSort = TRUE} then the returned matrix is sorted to be consistent with the factor loading matrix.
#' \item \strong{Rhat}: (Matrix) The (possibly sorted) reproduced correlation matrix.If \code{itemSort = TRUE} then the returned matrix is sorted to be consistent with the factor loading matrix.
#' \item \strong{Resid}: (Matrix) A (possibly sorted) residual matrix (R - Rhat) for the between battery correlations.
#' \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{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{loadings}: (Matrix) the factor loadings,
#' \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.
#' \item \strong{rotate} (Character) The chosen rotation algorithm.
#' \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{Call}: (call) A copy of the function call.
#' }
#'
#' @references Boruch, R. F., Larkin, J. D., Wolins, L., & MacKinney, A. C. (1970).
#' Alternative methods of analysis: Multitrait-multimethod data. \emph{Educational
#' and Psychological Measurement, 30}(4), 833–853.
#' https://doi.org/10.1177/0013164470030004055
#' @references Browne, M. W. (1979). The maximum-likelihood solution in inter-battery factor analysis.
#' \emph{British Journal of Mathematical and Statistical Psychology, 32}(1), 75-86.
#' @references Browne, M. W. (1980). Factor analysis of multiple batteries by maximum likelihood.
#' \emph{British Journal of Mathematical and Statistical Psychology, 33}(2), 184-199.
#' @references Browne, M. W. (2001). An overview of analytic rotation in
#' exploratory factor analysis. \emph{Multivariate Behavioral Research, 36}(1), 111-150.
#' @references Burnham, K. P. & Anderson, D. R. (2004). Multimodel inference: Understanding AIC and BIC in model selection.
#' \emph{Sociological methods and research, 33}, 261-304.
#' @references Cudeck, R. (1982). Methods for estimating between-battery factors,
#' \emph{Multivariate Behavioral Research, 17}(1), 47-68. 10.1207/s15327906mbr1701_3
#' @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 Tucker, L. R. (1958). An inter-battery method of factor analysis.
#' \emph{Psychometrika, 23}(2), 111-136.
#' @family Factor Analysis Routines
#'
#' @author
#' \itemize{
#' \item Niels G. Waller (nwaller@umn.edu)
#' \item Casey Giordano (Giord023@umn.edu)
#'}
#'
#' @examples
#'
#' # Example 1:
#' # Example from: Browne, M. W. (1979).
#' #
#' # Data originally reported in:
#' # Thurstone, L. L. & Thurstone, T. G. (1941). Factorial studies
#' # of intelligence. Psychometric Monograph (2), Chicago: Univ.
#' # Chicago Press.
#'
#' R.XY <- matrix(c(
#' 1.00, .554, .227, .189, .461, .506, .408, .280, .241,
#' .554, 1.00, .296, .219, .479, .530, .425, .311, .311,
#' .227, .296, 1.00, .769, .237, .243, .304, .718, .730,
#' .189, .219, .769, 1.00, .212, .226, .291, .681, .661,
#' .461, .479, .237, .212, 1.00, .520, .514, .313, .245,
#' .506, .530, .243, .226, .520, 1.00, .473, .348, .290,
#' .408, .425, .304, .291, .514, .473, 1.00, .374, .306,
#' .280, .311, .718, .681, .313, .348, .374, 1.00, .672,
#' .241, .311, .730, .661, .245, .290, .306, .672, 1.00), 9, 9)
#'
#'
#' dimnames(R.XY) <- list(c( paste0("X", 1:4),
#' paste0("Y", 1:5)),
#' c( paste0("X", 1:4),
#' paste0("Y", 1:5)))
#'
#' out <- faIB(R = R.XY,
#' n = 710,
#' NVarX = 4,
#' numFactors = 2,
#' itemSort = FALSE,
#' rotate = "oblimin",
#' rotateControl = list(standardize = "Kaiser",
#' numberStarts = 10),
#' Seed = 1)
#'
#' # Compare with Browne 1979 Table 2.
#' print(round(out$loadings, 2))
#' cat("\n\n")
#' print(round(out$Phi,2))
#' cat("\n\n MAD = ", round(out$fit["MAD"], 2),"\n\n")
#' print( round(out$facIndeterminacy,2) )
#'
#'
#' # Example 2:
#' ## Correlation values taken from Boruch et al.(1970) Table 2 (p. 838)
#' ## See also, Cudeck (1982) Table 1 (p. 59)
#' corValues <- c(
#' 1.0,
#' .11, 1.0,
#' .61, .47, 1.0,
#' .42, -.02, .18, 1.0,
#' .75, .33, .58, .44, 1.0,
#' .82, .01, .52, .33, .68, 1.0,
#' .77, .32, .64, .37, .80, .65, 1.0,
#' .15, -.02, .04, .08, .12, .11, .13, 1.0,
#' -.04, .22, .26, -.06, .07, -.10, .07, .09, 1.0,
#' .13, .21, .23, .05, .07, .06, .12, .64, .40, 1.0,
#' .01, .04, .01, .16, .05, .07, .05, .41, -.10, .29, 1.0,
#' .27, .13, .18, .17, .27, .27, .27, .68, .18, .47, .33, 1.0,
#' .24, .02, .12, .12, .16, .23, .18, .82, .08, .55, .35, .76, 1.0,
#' .20, .18, .16, .17, .22, .11, .29, .69, .20, .54, .34, .68, .68, 1.0)
#'
#' ## Generate empty correlation matrix
#' BoruchCorr <- matrix(0, nrow = 14, ncol = 14)
#'
#' ## Add upper-triangle correlations
#' BoruchCorr[upper.tri(BoruchCorr, diag = TRUE)] <- corValues
#' BoruchCorr <- BoruchCorr + t(BoruchCorr) - diag(14)
#'
#' ## Add variable names to the correlation matrix
#' varNames <- c("Consideration", "Structure", "Sup.Satisfaction",
#' "Job.Satisfaction", "Gen.Effectiveness", "Hum.Relations", "Leadership")
#'
#' ## Distinguish between rater X and rater Y
#' varNames <- paste0(c(rep("X.", 7), rep("Y.", 7)), varNames)
#'
#' ## Add row/col names to correlation matrix
#' dimnames(BoruchCorr) <- list(varNames, varNames)
#'
#' ## Estimate a model with one, two, and three factors
#' for (jFactors in 1:3) {
#' tempOutput <- faIB(R = BoruchCorr,
#' n = 111,
#' NVarX = 7,
#' numFactors = jFactors,
#' rotate = "oblimin",
#' rotateControl = list(standardize = "Kaiser",
#' numberStarts = 100))
#'
#' cat("\nNumber of inter-battery factors:", jFactors,"\n")
#' print( round(tempOutput$fit,2) )
#' } # END for (jFactors in 1:3)
#'
#' ## Compare output with Cudeck (1982) Table 2 (p. 60)
#' BoruchOutput <-
#' faIB(R = BoruchCorr,
#' n = 111,
#' NVarX = 7,
#' numFactors = 2,
#' rotate = "oblimin",
#' rotateControl = list(standardize = "Kaiser"))
#'
#' ## Print the inter-battery factor loadings
#' print(round(BoruchOutput$loadings, 3))
#' print(round(BoruchOutput$Phi, 3))
#'
#'
#'
#' @import GPArotation
#' @import stats
#' @export
faIB <-
function(X = NULL,
R = NULL,
n = NULL,
NVarX = 4,
numFactors = 2,
itemSort = FALSE,
rotate = "oblimin",
bootstrapSE = FALSE,
numBoot = 1000, ## Number of bootstrapped samples drawn
CILevel = .95, ## Bootstrap SE confidence level
rotateControl = NULL,
Seed = 1){
## Return a copy of the call function.
Call <- match.call()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# ----ERROR CHECKING ----
## 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))
## ____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) {
#----____END ERROR CHECKING----
## Set the seed for reproducibility
set.seed(Seed)
# Calculate R matrix if not given
if(!is.null(X)) R <- cor(X)
nVar <- ncol(R)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ##
#### DEFINE FUNCTIONS ####
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ##
## Function to compute model fit indices
ComputeModelFit <- function(R,
NVarX,
numFactors,
Lambda,
Phi){
## Define the number of variables in the model
nVar <- nrow(Lambda)
## Compute model-implied (reduced) correlation matrix
Rhat <- Lambda %*% Phi %*% t(Lambda)
## Find the residual correlation matrix
Resid <- R - Rhat
## Residual XY correlations
ResidRXY <- Resid[1:NVarX, (NVarX+1):nVar]
## Mean absolute difference of the residual XY correlation matrix
MAD <- mean( abs(ResidRXY) )
## Redefine the number of X variables
p <- NVarX
## Within X-battery correlation matrix
R.XX <- R[1:p, 1:p]
## Within Y-Battery correlation matrix
R.YY <- R[(p+1):nVar, (p+1):nVar]
## Between battery correlation matrix
R.XY <- R[1:p, (p+1):nVar]
# Squared canonical correlations
rhoSq <- Re(eigen(R.XY %*% solve(R.YY) %*% t(R.XY) %*% solve(R.XX))$values)
# Browne 1979 formula for corrected chi-square
chiSqTmp <- 1
for(i in (numFactors+1):NVarX ){
chiSqTmp <- chiSqTmp * (1 - rhoSq[i])
}
# see Browne 1979 p. 84
pstarInv <- 0
for(i in 1:numFactors){
pstarInv <- pstarInv + 1/rhoSq[i]
}
# Lawley 1959 correction for n
## NOTE: No correction if sample size (n) is not specified
if ( !is.null(n) ) {
nstar <- (n - 1 - numFactors - .5*(nVar+1) + pstarInv)
chiSq <- -nstar * log(chiSqTmp)
dfChiSq <- (NVarX - numFactors)*( (nVar - NVarX) - numFactors)
if(dfChiSq <= 0) warning("\n\n*****Zero Degrees of Freedom*****\n")
# check if rotation is orthogonal
# orthogonal is a logical (T/F)
orthogonal <- all.equal(Phi, diag(numFactors),
check.attributes = FALSE,
check.names = FALSE) == TRUE
# ~~~~~~~~~~~~~~AIC and BIC~~~~~~~~~~~~~~~~~
# qk = parsimony corrections factor
if (orthogonal) {
qk <- numEstParams <- (nVar * numFactors)
} else {
qk <- numEstParams <- (nVar * numFactors) + numFactors*(numFactors - 1)/2
} # END if (orthogonal)
#BIC
BIC <- chiSq + log(n) * qk
#AIC
AIC <- chiSq + 2 * qk
} else {
## If sample size is NULL, set the following as NA
BIC <- AIC <- chiSq <- dfChiSq <- NA
} # END if ( !is.null(n) )
## Return
list(Rhat = Rhat,
Resid = Resid,
MAD = MAD,
chiSq = chiSq,
DF = dfChiSq,
AIC = AIC,
BIC = BIC)
}# ED Compute Model Fit
## Extract names of the indicators to rename final output
GetVarNames <- function(){
varNames <- NULL
## If variable names are provided, retain them
if ( !is.null(X) ) varNames <- colnames(X)
if ( !is.null(R) ) varNames <- colnames(R)
if(is.null(varNames )) varNames <- paste0("V", 1:nVar)
varNames
}
# Main code to compute MLE multiple battery factor loadings
MBFact <- function(R, NVarX){
NVar <- ncol(R)
p <- NVarX
R.XX <- R[1:p, 1:p]
R.YY <- R[(p+1):NVar, (p+1):NVar]
R.XY <- R[1:p, (p+1):NVar]
F <- solve(R.XX) %*% R.XY %*% solve(R.YY) %*% t(R.XY)
VDV.F <- eigen(F)
V.F <- VDV.F$vectors
D.F <- VDV.F$values
G <- diag(sqrt(diag(t(V.F) %*% R.XX %*% V.F)))
B.X <- V.F %*% solve(G)
B.Y <- solve(R.YY)%*% t(R.XY) %*% B.X %*% diag(1/sqrt(D.F))
# Unrotated loadings for X (sqrt of sigular values)
LX <-(R.XX %*% B.X)[, 1:numFactors, drop = FALSE] %*%
diag(D.F^.25)[1:numFactors, 1:numFactors]
# Unrotated loadings for Y
LY <- R.YY %*% B.Y[, 1:numFactors, drop = FALSE] %*%
diag(D.F^.25)[1:numFactors, 1:numFactors]
# &&&& Why do we have to ask for the Real part of the complex number
urLoadings <- Re(rbind(LX, LY))
urLoadings
}#END MBFact
SortItems <- function(Lambda1, varNames, NVarX){
soutX <- faSort(Lambda1[1:NVarX, ],
phi = NULL,
BiFactor = FALSE,
salient = 0.25,
reflect = FALSE)
Lambda1[1:NVarX, ] <- soutX$loadings
soutY <- faSort(Lambda1[(NVarX + 1):nVar, ],
phi = NULL,
BiFactor = FALSE,
salient = 0.25,
reflect = FALSE)
Lambda1[(NVarX + 1):nVar, ] <- soutY$loadings
newOrder <- c(soutX$sortOrder,
(soutY$sortOrder + NVarX))
rownames(Lambda1) <- varNames[newOrder]
list( Lambda1 = Lambda1,
newOrder = newOrder)
}# END SortItems
#----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
## Compute Guttman's factor determinacy indices
GuttmanIndices <- function(Lambda,
PhiMat,
SampCorr) {
## Purpose: Compute Guttman (1955) factor indeterminacy indices
##
## Args: Lambda: (Matrix) Rotated factor loadings matrix
## PhiMat: (Matrix) Factor correlation matrix
## Fator structure (works for either oblique or orthogonal model)
facStruct <- Lambda %*% PhiMat
Rinv <-try(solve(SampCorr), silent = TRUE)
## If non-invertible, return NA values 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
## 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
##
## Args: lambda: (matrix) unrotated loadings to rotate
## rotation: (character) which rotation to use
## spinMatrix: (Matrix) randomly spin unrotated factor loadings
## 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" = {
list(loadings = lambda,
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" = {
promaxQ(urLoadings = lambda,
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
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# ---- MAIN PROGRAM----
# ----____Get varNames----
varNames <- GetVarNames()
#----____Compute MB factor loadings----
urLoadings <- MBFact(R, NVarX)
## If only 1 factor, rotate must equal "none"
if ( numFactors == 1 ) rotate <- "none"
## ----____ Standardize loadings ####
Stnd <- faStandardize(method = cnRotate$standardize,
lambda = urLoadings)
## Extract DvInv for later unstandardization
lambda <- Stnd$lambda
##----____ Rotate loadings from random spins ####
## 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(randSpinMatrix) {
rotatedLambda <- internalRotate(lambda = lambda,
rotation = rotate,
spinMatrix = randSpinMatrix,
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[]
## Return whole list, not just 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]]
# December 9, 2019 Note that order factors (Factor Sort) is not yet implimented
#### -------____ Sort rotated factor solutions ####
## 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]
## ----____Compute Factor Indeterminancy ####
## if the user provides urloadings then do not
## compute factor indeterminancy values
## Guttman's factor indeterminacy indices
uniqueSolutions[[iternumStart]]$facIndeterminacy <-
GuttmanIndices(Lambda = selected$loadings,
PhiMat = selected$Phi,
SampCorr = R ) # perhaps change name of R matrix
## 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 create solution sets
DisVal <- unlist(lapply(uniqueSolutions, function(x) x$RotationComplexityValue))
## Determine the number of digits for rotation convergence
## NOTE: subtract 2 because the leading zero and decimal count as characters
## NOTE: Convert scientific notation (e.g., 1e-5) to decimal form (works if non-scientific notation)
numDigits <- nchar( format(cnRotate$epsilon, scientific = FALSE) ) - 2
## Round the discrepancy values to specified number of digits
DisVal <- round(x = DisVal,
digits = numDigits)
## 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)
## Select the loading matrix with the lowest discrep value
## 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("FacIndeterminacy" = uniqueSolutions[[1]]$facIndeterminacy)
## If factor indeter. not computed, give NA values
if ( is.null(facIndeter) ) facIndeter <- rep(NA, ncol(minLambda))
Lambda1 <- minLambda
## ----____Reflect factors ----
# If reflect = TRUE then reflect factors s.t.
# salient loadings are positive
## Determine whether factors are negatively oriented
if (numFactors > 1) {
Dsgn <- diag(sign(colSums(Lambda1^3)))
} else {
Dsgn <- matrix(sign(colSums(Lambda1^3)), 1, 1)
} # END if (numFactors > 1)
## If factors negatively oriented, multiply by -1, else multiply by 1
Lambda1 <- Lambda1 %*% Dsgn
if (!is.null(minPhi)) {
## If factor is negative, reverse corresponding factor correlations
minPhi <- Dsgn %*% minPhi %*% Dsgn
} # END if (!is.null(minPhi))
# ----____Sort Items?----
# if no item sort, newOrder is original order
newOrder <- 1:nVar
if(itemSort == TRUE){
if(numFactors == 1){
newOrderX <- sort.list(Lambda1[1:NVarX])
newOrderY <- sort.list(Lambda1[(NVarX+1):nVar])
newOrder <- c(newOrderX , (newOrderY+NVarX) )
Lambda1 <- matrix(Lambda1[newOrder], nVar, 1)
rownames(Lambda1) <- varNames[newOrder]
colnames(Lambda1) <- "F1"
PhiNames <- colnames(Lambda1) <- paste0("F", 1:ncol(Lambda1))
dimnames(minPhi) <- list(PhiNames, PhiNames)
}
if(numFactors > 1){
sorted_out <- SortItems(Lambda1, varNames, NVarX)
Lambda1 <- sorted_out$Lambda1
newOrder <- sorted_out$newOrder
PhiNames <- colnames(Lambda1) <- paste0("F", 1:ncol(Lambda1))
dimnames(minPhi) <- list(PhiNames, PhiNames)
}
}else{
rownames(Lambda1) <- varNames
PhiNames <- colnames(Lambda1) <- paste0("F", 1:ncol(Lambda1))
dimnames(minPhi) <- list(PhiNames, PhiNames)
} #END Sort Items?
# if itemSort = TRUE the reOrder items in R
R <- R[newOrder, newOrder]
#---- ____Compute Model Fit----
fitOut <- ComputeModelFit(R = R,
NVarX = NVarX,
numFactors = numFactors,
Lambda = Lambda1,
Phi = minPhi)
Rhat <- fitOut$Rhat
Resid <- fitOut$Resid
MAD <- fitOut$MAD
chiSq <- fitOut$chiSq
DF <- fitOut$DF
AIC <- fitOut$AIC
BIC <- fitOut$BIC
#### -------- BOOTSTRAP SETUP -------- ####
## If true, compute bootstrap standard errors
if (bootstrapSE == TRUE) {
#This reorders the cols of X if itemSort = TRUE
X <- X[, newOrder]
# Number of subjects
nSubj <- nrow(X)
# Lists/variables to hold output
## Hold factor loadings from bootstraps
fList <-
## Hold factor correlations from bootstraps
phiList <-
## Hold factor indeterminacy indices from bootstraps
FIList <- 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, ])
## ____ Extract factor loadings ####
## Extract unrotated factors using resampled data matrix
bsLambda <- MBFact(R = Rsamp,
NVarX = NVarX)
## Conduct standardization
bsStnd <- faStandardize(method = cnRotate$standardize,
lambda = bsLambda)
## Extract the standardized bootstrapped (unrotated) factor loadings
bsLambda <- bsStnd$lambda
## ____ Rotate from random spins ####
## Find the "global" min out of all the random starting configurations
## Initialize a list to store rotations from random spins
bsStarts <- vector("list", cnRotate$numberStarts)
## Populate the list with random orthogonal matrices (i.e., start values for rotation)
bsStarts <- lapply(bsStarts, function(x) randStart(dimension = numFactors))
## Conduct rotations from random start values
bsStarts <- lapply(bsStarts, function(randSpinMatrix) {
## Rotate the bootstrapped samples
bsRotated <- internalRotate(lambda = bsLambda,
rotation = rotate,
spinMatrix = randSpinMatrix,
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
bsRotated
}) # END bsStarts <- lapply(bsStarts, function(x)
## ____ Order unique solutions ####
## 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
bsMinimum <- which.min(bootstrapComplexityFunc)
## Extract global minimum factor loadings and Phi matrices from
bsLambda <- bsStarts[[bsMinimum]]$loadings
bsPhi <- bsStarts[[bsMinimum]]$Phi
## Unstandardize the 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 = Lambda1,
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^3) < 0 ) bsLambda <- bsLambda * -1
## Define newly aligned elements
AlignedLambda <- bsLambda
AlignedPhi <- bsPhi
} # END if (numFactors == 1)
## Save loadings as the bootstrapped factor loadings
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)
} # END for (iSample in seq_len(numBoot))
#### -----____ BootStrap stnd 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))
## Bootstrap standard errors for factor loadings
loadSE <- apply(loadArray, 1:2, sd)
# Bootstrap standard errors for factor correlations
phiSE <- apply(phiArray, 1:2, sd)
## ____ Confidence intervals ####
## Set alpha level
alpha <- 1 - CILevel
## Confidence interval matrices for factor loadings
loadCI.upper <- apply(loadArray, 1:2, quantile, 1 - (alpha / 2))
loadCI.lower <- apply(loadArray, 1:2, quantile, (alpha / 2))
# 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))
## ____ Bootstrap Fac Indeterminacy SE ####
## If none of the fac indeter values are NA, proceed. Else, set as NA
if (!any(is.na(FIArray))) {
## Compute the standard error of the bootstrap factor indetermancy values
FISE <- apply(FIArray, 1:2, sd)
## Find the upper/lower bounds of the confidence interval
FICI.upper <- apply(FIArray, 2, quantile, 1 - (alpha / 2))
FICI.lower <- apply(FIArray, 2, quantile, (alpha / 2))
## Create data frame of fac indeterminacy estimates
facIndeter <- data.frame(facIndeter,
t(FISE),
FICI.lower,
FICI.upper)
## 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"))
## Add row names to correspond to each factor
rownames(facIndeter) <- paste0("f", 1:numFactors)
}else{
facIndeter <- NA
FISE <- NA
} # END if (!any(is.na(FIArray)))
#### ------- NAME BOOTSTRAP OBJECTS -------- ####
## Dimension names
## Name the factor indicators
rownames(loadSE) <-
rownames(loadCI.upper) <-
rownames(loadCI.lower) <- varNames[newOrder]
## Name the factors
colnames(loadSE) <-
colnames(loadCI.upper) <-
colnames(loadCI.lower) <-
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) {
loadSE <-
loadCI.upper <-
loadCI.lower <-
loadArray <-
phiSE <-
phiCI.upper <-
phiCI.lower <-
phiArray <-
FIArray <-
FISE <-
FICI.upper <-
FICI.lower <- NULL
} # END if (bootstrapSE == FALSE)
## CG Edits (13 dec 2019): If sample size is NULL, pchisq cannot be computed
if ( !is.null(n) ) {
fitVec <- c(chiSq, DF, pchisq(chiSq, DF), MAD, AIC, BIC)
} else {
fitVec <- c(chiSq, DF, NA, MAD, AIC, BIC)
} # END if ( !is.null(n) )
names(fitVec) <- c("chiSq", "DF", "p-value", "MAD", "AIC", "BIC")
#----____Output ----
invisible(list(loadings = Lambda1, # this is min complexity solutution
Phi = minPhi,
fit = fitVec,
R = R,
Rhat = Rhat,
Resid = Resid,
facIndeterminacy = facIndeter,
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,
rotate = rotate,
loadingsArray = loadArray,
PhiArray = phiArray,
facIndeterminacyArray = FIArray,
rotateControl = cnRotate,
unSpunSolution = uniqueSolutions[[UnSpunPosition]],
Call = Call
))#END inivisble(list( . . .
} ##END faMB
Try the fungible package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.