Nothing
#' Regularized Factor Analysis
#'
#' This function applies the regularized factoring method to extract an unrotated factor structure matrix.
#'
#' @param R (Matrix) A correlation matrix to be analyzed.
#' @param numFactors (Integer) The number of factors to extract. Default: numFactors = 1.
#' @param facMethod (Character) "rls" for regularized least squares estimation or
#' "rml" for regularized maximum likelihood estimation. Default: facMethod = "rls".
#' @return The main output is the matrix of unrotated factor loadings.
#' \itemize{
#' \item \strong{loadings}: (Matrix) A matrix of unrotated factor loadings.
#' \item \strong{h2}: (Vector) A vector of estimated communality values.
#' \item \strong{L}: (Numeric) Value of the estimated penality parameter.
#' \item \strong{Heywood} (Logical) TRUE if a Heywood case is detected
#' (this should never happen).
#' }
#'
#' @author
#' Niels G. Waller (nwaller@umn.edu)
#'
#' @family Factor Analysis Routines
#'
#' @references Jung, S. & Takane, Y. (2008). Regularized common factor analysis.
#' New trends in psychometrics, 141-149.
#'
#' @examples
#'
#' data("HW")
#'
#' # load first HW data set
#'
#' RHW <- cor(x = HW$HW6)
#'
#' # Compute principal axis factor analysis
#' fapaOut <- faMain(R = RHW,
#' numFactors = 3,
#' facMethod = "fapa",
#' rotate = "oblimin",
#' faControl = list(treatHeywood = FALSE))
#'
#'
#' fapaOut$faFit$Heywood
#' round(fapaOut$h2, 2)
#'
#' # Conduct a regularized factor analysis
#' regOut <- fareg(R = RHW,
#' numFactors = 3,
#' facMethod = "rls")
#' regOut$L
#' regOut$Heywood
#'
#'
#' # rotate regularized loadings and align with
#' # population structure
#' regOutRot <- faMain(urLoadings = regOut$loadings,
#' rotate = "oblimin")
#'
#' # ALign
#' FHW <- faAlign(HW$popLoadings, fapaOut$loadings)$F2
#' Freg <- faAlign(HW$popLoadings, regOutRot$loadings)$F2
#'
#' AllSolutions <- round(cbind(HW$popLoadings, Freg, FHW),2)
#' colnames(AllSolutions) <- c("F1", "F2", "F3", "Fr1", "Fr2", "Fr3",
#' "Fhw1", "Fhw2", "Fhw3")
#' AllSolutions
#'
#'
#' rmsdHW <- rmsd(HW$popLoadings, FHW,
#' IncludeDiag = FALSE,
#' Symmetric = FALSE)
#'
#' rmsdReg <- rmsd(HW$popLoadings, Freg,
#' IncludeDiag = FALSE,
#' Symmetric = FALSE)
#'
#' cat("\nrmsd HW = ", round(rmsdHW,3),
#' "\nrmsd reg = ", round(rmsdReg,3))
#'
#'
#'
#' @export
#--------------------------------------------------------#
fareg <- function(R,
numFactors = 1,
facMethod = "rls"){
# Date: 04/19/2019
# Arguments:
#
# R (numeric matrix) Input correlation matrix.
# numfactors (scaler) Number of factors to extract, Default = 1
# facMethod (Character) "rls" for regulaized least squares estimation or
# "rml" for regularixed maximum likelihood estimation.
#Values
# loadings,
# h2,
# L,
# Heywood
if(!isSymmetric(R)) stop("\nInput data must be a correlation matrix")
if( !all.equal( as.vector(diag(R)), rep(1, dim(R)[1]) ) ){
stop("\nInput data must be a correlation matrix")
}
PlausibleMethods <- c("rls", "rml")
if(facMethod %in% PlausibleMethods == FALSE){
stop("\nInvalid argument for facMethod")
}
p <- nrow(R)
communality <- "SMC"
## Check to ensure that R is positive definite
eigenVal <- eigen(R)$value
if ( min(eigenVal) <= 1E-6 ) {
warning("Inverting the correlation matrix for SMC communality estimates requires a positive-definite matrix.")
communality <- "maxr"
} # END if ( min(eigenVal) <= 1E-8 )
if(communality == "SMC"){
Psi <- diag(diag(solve(R))^-1)
}
if(communality == "maxr"){
R0 <- R
diag(R0) <- 0
h2temp <- apply(abs(R0),2, max)
Psi <- diag(1 - h2temp)
}
#regularization parameter
Lstart = 1
# function to minimize
FncRidgeULS <- function(L){
Ri <- R - (L * Psi)
UDU <- eigen(Ri)
U <- UDU$vectors
D <- UDU$values
if( numFactors == 1){
F <- U[,1] * sqrt(D[1])
}
if( numFactors > 1){
F <- U[,1:numFactors] %*% diag(sqrt(D[1:numFactors]))
}
Rhat <- F %*% t(F)
#return fnc
rmsd(Ri, Rhat,
Symmetric = TRUE,
IncludeDiag = FALSE)
}
# function to minimize
FncRidgeMLE <- function(L){
Ri <- R - (L * Psi)
UDU <- eigen(Ri)
U <- UDU$vectors
D <- UDU$values
if( numFactors == 1){
F <- U[,1] * sqrt(D[1])
}
if( numFactors > 1){
F <- U[,1:numFactors] %*% diag(sqrt(D[1:numFactors]))
}
Rhat <- F %*% t(F)
diag(Rhat) <- 1
RhatInv <- solve(Rhat)
log(det(Rhat)) + sum(diag(R %*% RhatInv )) - log(det(R)) - p
}
# function to produce loadings from optimal L
faSolution <- function(L){
Ri <- R - (L * Psi)
UDU <- eigen(Ri)
U <- UDU$vectors
D <- UDU$values
# return unrotated factor loading matrix
if( numFactors == 1){
F <- U[,1] * sqrt(D[1])
}
if( numFactors > 1){
F <- U[,1:numFactors] %*% diag(sqrt(D[1:numFactors]))
}
F
}
maxPsi <- max(diag(Psi))
Lup <- 1/maxPsi + .001
# optimize function
if(facMethod == "rls"){
out <- optimize(f = FncRidgeULS, lower = .001, upper = Lup)
}
if(facMethod == "rml"){
out <- optimize(f = FncRidgeMLE, lower = 0, upper = Lup)
}
Lopt <- out$minimum
# set convergence status
converged <- FALSE
if(Lopt > 0 & Lopt < Lup) converged <- TRUE
# Compute regularized loadings
loadings <- faSolution(Lopt)
if(numFactors == 1){
loadings <- as.matrix(loadings)
h2 <- as.vector(loadings^2)
}
if( numFactors > 1){
h2 <- apply(loadings^2, 1, sum)
}
Heywood <- FALSE
# By design, Heywood should never be TRUE
if(max(h2 > 1)) Heywood <- TRUE
list(loadings = loadings,
h2 = h2,
L = Lopt,
converged = converged,
Heywood = Heywood)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.