# --- Basic rotation calls ---
data(Harman, package = "GPArotation") # 8 physical variables
quartimax(Harman8) # direct rotation call
GPFRSorth(Harman8, method = "quartimax") # equivalent via wrapper
GPFRSoblq(Harman8, method = "quartimin", normalize = TRUE)
loadings(quartimin(Harman8, normalize = TRUE)) # extract loadings directly
# --- Passing criterion arguments via methodArgs ---
# Crawford-Ferguson family: kappa selects the criterion.
# For box26: p = 26 variables, m = 3 factors.
# Equamax: kappa = m / (2 * p) = 3 / 52
# Parsimax: kappa = (m - 1) / (p + m - 2) = 2 / 27
data(Thurstone, package = "GPArotation") # 26 variable box problem
GPFRSoblq(box26, method = "cf", methodArgs = list(kappa = 3/52)) # Equamax
GPFRSoblq(box26, method = "cf", methodArgs = list(kappa = 2/27)) # Parsimax
# --- Two-step vs single-step factanal for oblique rotation ---
#
# The recommended approach for oblique rotation is the two-step procedure:
# (1) obtain unrotated loadings from factanal, then
# (2) rotate separately using GPArotation.
# This gives full control over the rotation, including random starts.
#
# Prior to R 4.5.1, the single-step approach (rotation inside factanal)
# had a bug in factor reordering after oblique rotation.
# This was fixed by the R core team in R 4.5.1.
data("WansbeekMeijer", package = "GPArotation")
# Step 1: unrotated 3-factor solution
fa.unrotated <- factanal(factors = 3, covmat = NetherlandsTV,
normalize = TRUE, rotation = "none")
# Step 2: oblique Crawford-Ferguson rotation with kappa = 0.3
# (non-standard kappa, not corresponding to any named special case)
set.seed(44)
fa.cf <- cfQ(fa.unrotated, kappa = 0.3, normalize = TRUE, randomStarts = 100)
fa.cf
# Single-step via factanal - correct in R >= 4.5.1
if (getRversion() >= "4.5.1") {
set.seed(44)
fa.factanal <- factanal(factors = 3, covmat = NetherlandsTV, rotation = "cfQ",
control = list(rotate = list(normalize = TRUE, kappa = 0.3, randomStarts = 100)))
# The two approaches should agree after sorting
fa.sorted <- print(fa.cf, sortLoadings = TRUE)
cat("Maximum difference in loadings between two-step and single-step:\n")
print(max(abs(abs(fa.sorted$loadings) - abs(fa.factanal$loadings))))
} else {
cat("Single-step factanal oblique rotation requires R >= 4.5.1.\n")
cat("Use the two-step procedure above for correct results.\n")
}
# --- Displaying rotation output ---
origdigits <- options("digits")
data("CCAI", package = "GPArotation")
fa.unrotated <- factanal(factors = 3, covmat = CCAI_R, n.obs = 461, rotation = "none")
res <- oblimin(fa.unrotated, gam = -0.5, randomStarts = 20)
# gam = -0.5: more orthogonal than quartimin
res # default print
print(res) # equivalent to above
print(res, Table = TRUE) # include iteration table
print(res, rotateMat = TRUE) # include rotating matrix
print(res, digits = 2) # rounded to 2 decimal places
summary(res) # pattern and structure matrices for oblique rotation
summary(res, Structure = FALSE) # pattern matrix only
options(digits = origdigits$digits)
# \donttest{
# --- Random start diagnostics ---
# When randomStarts > 1, the output includes randStartChar which summarizes
# the random start results:
# randomStarts : number of random starts attempted
# Converged : number of starts that converged
# atMinimum : number of starts at the same lowest minimum
# localMins : number of distinct local minima found
data(Thurstone, package = "GPArotation")
res <- GPFRSoblq(box26, method = "geomin", normalize = TRUE, randomStarts = 50)
res$randStartChar
# --- Factor ordering ---
# Raw GPArotation output is unsorted — factors may appear in any order
# depending on the starting matrix. Use print() to obtain sorted loadings.
# Once sorted, repeated calls to print() are stable.
set.seed(334)
xusl <- quartimin(Harman8, normalize = TRUE, randomStarts = 100)
loadings(xusl) # unsorted raw output
max(abs(print(xusl)$loadings - xusl$loadings)) == 0 # FALSE: print() reorders
xsl <- print(xusl) # capture sorted result
max(abs(print(xsl)$loadings - xsl$loadings)) == 0 # TRUE: already sorted
# --- Normalization ---
# Kaiser normalization
data("CCAI", package = "GPArotation")
fa.unrotated <- factanal(factors = 3, covmat = CCAI_R, n.obs = 461, rotation = "none")
oblimin(fa.unrotated, normalize = TRUE, randomStarts = 100)
data("CCAI", package = "GPArotation")
res.u <- factanal(covmat = CCAI_R, factors = 4, rotation = "none",
n.obs = 461)
# TandemI --- non-smooth orthogonal, BB and Cayley both faster
res.t1l <- tandemI(res.u, algorithm = "legacy", fwindow = 1, maxit = 10000)
res.t1b <- tandemI(res.u, algorithm = "bb", fwindow = 10)
res.t1c <- tandemI(res.u, algorithm = "cayley", fwindow = 10)
cat("TandemI legacy:", nrow(res.t1l$Table) - 1, "iterations\n")
cat("TandemI bb: ", nrow(res.t1b$Table) - 1, "iterations",
" max diff:", format(max(abs(res.t1l$loadings - res.t1b$loadings)),
scientific = TRUE, digits = 2), "\n")
cat("TandemI cayley:", nrow(res.t1c$Table) - 1, "iterations",
" max diff:", format(max(abs(res.t1l$loadings - res.t1c$loadings)),
scientific = TRUE, digits = 2), "\n")
# Entropy --- orthogonal, BB and Cayley both faster
res.t2l <- entropy(res.u, algorithm = "legacy", fwindow = 1, maxit = 10000)
res.t2b <- entropy(res.u, algorithm = "bb", fwindow = 10)
res.t2c <- entropy(res.u, algorithm = "cayley", fwindow = 10)
cat("Entropy legacy:", nrow(res.t2l$Table) - 1, "iterations\n")
cat("Entropy bb: ", nrow(res.t2b$Table) - 1, "iterations",
" max diff:", format(max(abs(res.t2l$loadings - res.t2b$loadings)),
scientific = TRUE, digits = 2), "\n")
cat("Entropy cayley:", nrow(res.t2c$Table) - 1, "iterations",
" max diff:", format(max(abs(res.t2l$loadings - res.t2c$loadings)),
scientific = TRUE, digits = 2), "\n")
# Oblimin --- smooth oblique, BB faster, same solution
res.ol <- oblimin(res.u, algorithm = "legacy", fwindow = 1, maxit = 10000)
res.ob <- oblimin(res.u, algorithm = "bb", fwindow = 10)
cat("Oblimin legacy:", nrow(res.ol$Table) - 1, "iterations\n")
cat("Oblimin bb: ", nrow(res.ob$Table) - 1, "iterations",
" max diff:", format(max(abs(res.ol$loadings - res.ob$loadings)),
scientific = TRUE, digits = 2), "\n")
# Simplimax --- non-smooth oblique with many local minima.
# BB finds a substantially better solution in fewer iterations.
# Difference in loadings reflects different local minima, not noise.
res.sl <- simplimax(res.u, algorithm = "legacy", fwindow = 1, maxit = 10000)
res.sb <- simplimax(res.u, algorithm = "bb", fwindow = 10)
cat("Simplimax legacy:", nrow(res.sl$Table) - 1, "iterations",
" f =", round(res.sl$Table[nrow(res.sl$Table), 2], 6), "\n")
cat("Simplimax bb: ", nrow(res.sb$Table) - 1, "iterations",
" f =", round(res.sb$Table[nrow(res.sb$Table), 2], 6), "\n")
cat("BB found a lower criterion value --- better solution quality.\n")
# }
Run the code above in your browser using DataLab