GPA1guide.R
In GPArotation: Gradient Projection Factor Rotation

### R code from vignette source 'GPA1guide.Stex'

###################################################
### code chunk number 1: GPA1guide.Stex:22-24
###################################################
 options(continue="  ")
pdf.options(pointsize = 8)


###################################################
### code chunk number 2: GPA1guide.Stex:44-45
###################################################
library("GPArotation")


###################################################
### code chunk number 3: GPA1guide.Stex:89-100
###################################################
data("Harman", package = "GPArotation")

# Calling a rotation directly
qHarman <- quartimax(Harman8)

# Equivalently, via the wrapper function
qHarman <- GPFRSorth(Harman8, method = "quartimax")

# Two equivalent ways to access the rotated loadings
loadings(qHarman)       # via extractor function (recommended)
qHarman$loadings        # via direct list access


###################################################
### code chunk number 4: GPA1guide.Stex:112-120
###################################################
data("CCAI", package = "GPArotation")
y <- factanal(factors = 3, covmat = CCAI_R, n.obs = 461, rotation = "none")
y.quart <- quartimax(y$loadings)
max(loadings(y.quart) %*% t(y.quart$Th) - loadings(y))
y.obli <- oblimin(y$loadings, normalize = TRUE, randomStarts = 15)
max(loadings(y.obli) %*% t(y.obli$Th) - loadings(y))
# last equation on Page 678
max(loadings(y.obli) - loadings(y) %*% solve(t(y.obli$Th)))


###################################################
### code chunk number 5: GPA1guide.Stex:126-129
###################################################
y <- factanal(factors = 3, covmat = CCAI_R, n.obs = 461, rotation = "none")
y.obli <- oblimin(y$loadings, normalize = TRUE, randomStarts = 15)
max(abs(y.obli$Phi - t(y.obli$Th) %*% y.obli$Th))


###################################################
### code chunk number 6: GPA1guide.Stex:143-154
###################################################
set.seed(334)
res <- quartimin(Harman8, normalize = TRUE, randomStarts = 100)

# Raw unsorted loadings
loadings(res)

# Sorted loadings via print (factors reordered and signs adjusted)
res.sorted <- print(res)

# Once sorted, repeated calls to print are stable
max(abs(print(res.sorted)$loadings - res.sorted$loadings)) == 0  # TRUE


###################################################
### code chunk number 7: GPA1guide.Stex:193-198
###################################################
res.obli <- oblimin(Harman8, normalize = TRUE, randomStarts = 100)
# Pattern matrix (unsorted)
summary(res.obli, Structure = FALSE)
# Structure matrix
summary(res.obli, Structure = TRUE)


###################################################
### code chunk number 8: GPA1guide.Stex:241-261
###################################################
plotPatternStructure <- function(pattern, structure,
                                  labels = NULL,
                                  main = "Pattern vs Structure") {
  k   <- ncol(pattern)
  col <- palette.colors(k, palette = "Okabe-Ito")
  par(mfrow = c(1, k))
  for (j in 1:k) {
    lims <- range(c(pattern[, j], structure[, j]))
    plot(pattern[, j], structure[, j],
         xlim = lims, ylim = lims,
         xlab = "Pattern loading",
         ylab = "Structure loading",
         main = paste(if (!is.null(labels)) labels[j]
                      else paste("Factor", j)),
         pch  = 19, col = col[j])
    abline(0, 1, lty = 2, col = "grey60")
    abline(h = 0, col = "grey80")
    abline(v = 0, col = "grey80")
  }
}


###################################################
### code chunk number 9: GPA1guide.Stex:264-270
###################################################
res.obli.s  <- print(res.obli)
Pattern     <- loadings(res.obli.s)
Structure   <- Pattern %*% res.obli.s$Phi

plotPatternStructure(Pattern, Structure,
                     labels = c("Factor 1", "Factor 2"))


###################################################
### code chunk number 10: GPA1guide.Stex:295-299
###################################################
data(Thurstone, package = "GPArotation")
infomaxQ(box26, randomStarts = 100)         # 100 random starts
infomaxQ(box26, Tmat = Random.Start(3))     # single random start
infomaxQ(box26, randomStarts = 1)           # also a single random start


###################################################
### code chunk number 11: GPA1guide.Stex:328-332
###################################################
res <- geominQ(box26, normalize = TRUE, randomStarts = 100)
res$randStartChar
# Criterion value at best solution
res$Table[nrow(res$Table), "f"]


###################################################
### code chunk number 12: GPA1guide.Stex:357-413
###################################################
plotRotationLandscape <- function(A, method = "quartimax", n = 1000,
                                   main = NULL, ...) {
  # Plot the objective function landscape for 2-factor orthogonal rotation.
  # For 2 factors, all orthogonal rotations are parameterized by a single
  # angle theta in [0, 2*pi), giving a clean 1D landscape.
  #
  # Args:
  #   A      : a 2-factor unrotated loading matrix
  #   method : rotation criterion (default "quartimax")
  #   n      : number of angles to evaluate (default 1000)
  #   main   : plot title (default: "Rotation landscape: <method>")
  #   ...    : additional arguments passed to the vgQ criterion function

  if (ncol(A) != 2)
    stop("plotRotationLandscape only works for 2-factor solutions.")

  vgQfun_fn <- get(paste("vgQ", method, sep = "."),
                   envir = asNamespace("GPArotation"))

  if (is.null(main))
    main <- paste("Rotation landscape:", method)

  theta  <- seq(0, 2 * pi, length.out = n)
  f_vals <- numeric(n)

  for (i in seq_along(theta)) {
    Tmat      <- matrix(c( cos(theta[i]), sin(theta[i]),
                          -sin(theta[i]), cos(theta[i])), 2, 2)
    L         <- A %*% Tmat
    VgQ       <- do.call(vgQfun_fn, append(list(L), list(...)))
    f_vals[i] <- VgQ$f
  }

  # Find global minimum
  min_idx <- which.min(f_vals)

  plot(theta, f_vals,
       type  = "l",
       lwd   = 2,
       main  = main,
       xlab  = expression(theta ~ "(radians)"),
       ylab  = "f",
       xaxt  = "n")
  axis(1, at     = c(0, pi/2, pi, 3*pi/2, 2*pi),
           labels = c("0", expression(pi/2), expression(pi),
                      expression(3*pi/2), expression(2*pi)))
  abline(h   = f_vals[min_idx], col = "grey80", lty = 2)
  points(theta[min_idx], f_vals[min_idx],
         col = "tomato", pch = 19, cex = 1.5)
  text(theta[min_idx], f_vals[min_idx],
       labels = paste0("min f = ", round(f_vals[min_idx], 4)),
       pos    = 4,
       cex    = 0.8)

  invisible(data.frame(theta = theta, f = f_vals))
}


###################################################
### code chunk number 13: GPA1guide.Stex:419-420
###################################################
data(Harman, package = "GPArotation")


###################################################
### code chunk number 14: GPA1guide.Stex:423-428
###################################################
par(mfrow = c(2, 2))
plotRotationLandscape(Harman8, method = "quartimax")
plotRotationLandscape(Harman8, method = "varimax")
plotRotationLandscape(Harman8, method = "bentler")
plotRotationLandscape(Harman8, method = "entropy")


###################################################
### code chunk number 15: GPA1guide.Stex:455-468
###################################################
res <- plotRotationLandscape(Harman8, method = "simplimax", k = 4)

# Find all local minima by sign changes in the discrete derivative
df         <- diff(res$f)
local_mins <- which(df[-length(df)] < 0 & df[-1] > 0)

# Mark all local minima on the existing plot
points(res$theta[local_mins], res$f[local_mins],
       col = "steelblue", pch = 19, cex = 0.8)
legend("topright",
       legend = c("global minimum", "local minima"),
       col    = c("tomato", "steelblue"),
       pch    = 19, bty = "n", cex = 0.8)


###################################################
### code chunk number 16: GPA1guide.Stex:471-474
###################################################
cat("Total minima found:              ", length(local_mins), "\n")
cat("Expected due to symmetry:        ", 4, "\n")
cat("Approximate genuine local minima:", max(0, length(local_mins) - 4), "\n")


###################################################
### code chunk number 17: GPA1guide.Stex:505-536
###################################################
data(Harman, package = "GPArotation")
res.quart   <- quartimax(Harman8)
res.oblimin <- oblimin(Harman8)

L.quart   <- abs(loadings(res.quart))
L.oblimin <- abs(loadings(res.oblimin))

ord.quart   <- order(L.quart[, 1],   decreasing = TRUE)
ord.oblimin <- order(L.oblimin[, 1], decreasing = TRUE)

par(mfrow = c(1, 2), mar = c(5, 4, 4, 2))

barplot(t(L.quart[ord.quart, ]),
        beside      = TRUE,
        ylim        = c(0, 1),
        main        = "Quartimax",
        ylab        = "Absolute loading",
        xlab        = "Variable (sorted by Factor 1)",
        legend.text = c("Factor 1", "Factor 2"),
        args.legend = list(x = "topright"),
        col         = c("steelblue", "tomato"))

barplot(t(L.oblimin[ord.oblimin, ]),
        beside      = TRUE,
        ylim        = c(0, 1),
        main        = "Oblimin",
        ylab        = "Absolute loading",
        xlab        = "Variable (sorted by Factor 1)",
        legend.text = c("Factor 1", "Factor 2"),
        args.legend = list(x = "topright"),
        col         = c("steelblue", "tomato"))


###################################################
### code chunk number 18: GPA1guide.Stex:556-621
###################################################
plotSortedLoadings <- function(..., labels = NULL, col = NULL, 
                                main = "Sorted Absolute Loadings",
                                ylab = "Absolute loading", 
                                xlab = "Rank") {
  # Plot sorted absolute loadings for one or more GPArotation objects.
  # Multiple solutions are overlaid on a single plot for comparison.
  # Loadings are sorted from smallest to largest (left to right).
  #
  # Args:
  #   ...    : one or more GPArotation objects
  #   labels : character vector of legend labels (default: "Solution 1", etc.)
  #   col    : character vector of colors (default: auto-assigned)
  #   main   : plot title
  #   ylab   : y-axis label
  #   xlab   : x-axis label

  solutions <- list(...)

  for (i in seq_along(solutions)) {
    if (!inherits(solutions[[i]], "GPArotation"))
      stop("Argument ", i, " is not a GPArotation object.")
  }

  n <- length(solutions)

  if (is.null(labels))
    labels <- paste("Solution", seq_len(n))
  if (is.null(col))
    col <- palette.colors(n, palette = "Okabe-Ito")

  sorted_loadings <- lapply(solutions, function(x)
    sort(abs(as.vector(x$loadings)), decreasing = FALSE))

  all_values <- unlist(sorted_loadings)
  max_len    <- max(sapply(sorted_loadings, length))

  plot(NULL,
       xlim = c(1, max_len),
       ylim = c(0, max(all_values)),
       main = main,
       xlab = xlab,
       ylab = ylab,
       las  = 1)

  abline(h = seq(0, 1, by = 0.1), col = "grey90", lty = 1)

  for (i in seq_len(n)) {
    lines(seq_along(sorted_loadings[[i]]), sorted_loadings[[i]],
          col = col[i], lwd = 2)
    points(seq_along(sorted_loadings[[i]]), sorted_loadings[[i]],
           col = col[i], pch = 19, cex = 0.6)
  }

  legend("topleft", legend = labels, col = col, lwd = 2, pch = 19,
         bty = "n")

  invisible(sorted_loadings)
}

# Example
data(Harman, package = "GPArotation")
res.quart   <- quartimax(Harman8)
res.oblimin <- oblimin(Harman8)
res.geomin  <- geominT(Harman8)



###################################################
### code chunk number 19: GPA1guide.Stex:626-628
###################################################
plotSortedLoadings(res.quart, res.oblimin, res.geomin,
                   labels = c("Quartimax", "Oblimin", "Geomin"))


###################################################
### code chunk number 20: GPA1guide.Stex:685-716
###################################################
origdigits <- options("digits")
options(digits = 2)
trBritain <- matrix(c(.783,-.163,.811,.202,.724,.209,.850,.064,
  -.031,.592,-.028,.723,.388,.434,.141,.808,.215,.709),
  byrow = TRUE, ncol = 2)
trGermany <- matrix(c(.778,-.066,.875,.081,.751,.079,.739,.092,
  .195,.574,-.030,.807,-.135,.717,.125,.738,.060,.691),
  byrow = TRUE, ncol = 2)
# orthogonal rotation of trGermany towards trBritain
trx <- targetT(trGermany, Target = trBritain)
# Factor loadings after target rotation
trx
# Differences between loadings matrices after rotation
y <- trx$loadings - trBritain
print(y, digits = 1)
# Square root of the mean squared difference per item
sqrt(apply((y^2), 1, mean))
# Square root of the mean squared difference per factor
sqrt(apply((y^2), 2, mean))
# Identity coefficient per factor after rotation
2 * colSums(trx$loadings * trBritain) /
  (colSums(trx$loadings^2) + colSums(trBritain^2))
# Additivity coefficient per factor after rotation
diag(2 * cov(trx$loadings, trBritain)) /
  diag(var(trx$loadings) + var(trBritain))
# Proportionality coefficient per factor after rotation
colSums(trBritain * trx$loadings) /
  sqrt(colSums(trBritain^2) * colSums(trx$loadings^2))
# Correlation for each factor after rotation
diag(cor(trBritain, trx$loadings))
options(digits = origdigits$digits)


###################################################
### code chunk number 21: GPA1guide.Stex:728-751
###################################################
plot(trBritain[, 1], trBritain[, 2],
     xlim = c(-0.3, 1.0), ylim = c(-0.3, 1.0),
     xlab = "Factor 1", ylab = "Factor 2",
     main = "Target Rotation: Germany towards Britain",
     pch = 19, col = "steelblue", cex = 1.2)
abline(h = 0, lty = 2, col = "grey70")
abline(v = 0, lty = 2, col = "grey70")
points(trGermany[, 1], trGermany[, 2],
       pch = 17, col = "tomato", cex = 1.2)
points(loadings(trx)[, 1], loadings(trx)[, 2],
       pch = 15, col = "orange", cex = 1.2)
for (i in 1:nrow(trGermany)) {
  arrows(trGermany[i, 1], trGermany[i, 2],
         loadings(trx)[i, 1], loadings(trx)[i, 2],
         length = 0.08, col = "grey60")
}
legend("topright",
       legend = c("Britain (varimax rotated)", 
       		"East Germany (varimax rotated)", 
		"East Germany (rotated towards Britain)"),
       col    = c("steelblue", "tomato", "orange"),
       pch    = c(19, 17, 15),
       bty    = "n")


###################################################
### code chunk number 22: GPA1guide.Stex:778-792
###################################################
A <- matrix(c(.664, .688, .492, .837, .705, .82, .661, .457, .765, .322,
  .248, .304, -0.291, -0.314, -0.377, .397, .294, .428, -0.075, .192, .224,
  .037, .155, -.104, .077, -.488, .009), ncol = 3)
# using targetT
SPA <- matrix(c(rep(NA, 6), .7, .0, .7, rep(0, 3), rep(NA, 7),
  0, 0, NA, 0, rep(NA, 4)), ncol = 3)
xt <- targetT(A, Target = SPA)
# using pstT
SPApst <- matrix(c(rep(0, 6), .7, .0, .7, rep(0, 3), rep(0, 7),
  0, 0, 0, 0, rep(0, 4)), ncol = 3)
SPAW <- matrix(c(rep(0, 6), rep(1, 6), rep(0, 7), 1, 1, 0, 1,
  rep(0, 4)), ncol = 3)
xpst <- pstT(A, Target = SPApst, W = SPAW)
max(abs(loadings(xt) - loadings(xpst)))


###################################################
### code chunk number 23: GPA1guide.Stex:810-814
###################################################
data("CCAI", package = "GPArotation")
factanal(factors = 3, covmat = CCAI_R, n.obs = 461, rotation = "infomaxT")
factanal(factors = 3, covmat = CCAI_R, n.obs = 461, rotation = "infomaxT",
  control = list(rotate = list(normalize = TRUE, eps = 1e-6)))


###################################################
### code chunk number 24: GPA1guide.Stex:822-826
###################################################
data("WansbeekMeijer", package = "GPArotation")
fa.unrotated <- factanal(factors = 3, covmat = NetherlandsTV,
                         normalize = TRUE, rotation = "none")
quartimin(loadings(fa.unrotated), normalize = TRUE)


###################################################
### code chunk number 25: GPA1guide.Stex:836-859
###################################################
data("WansbeekMeijer", package = "GPArotation")
fa.unrotated <- factanal(factors = 3, covmat = NetherlandsTV,
                         normalize = TRUE, rotation = "none")
# Two-step procedure (always correct)
set.seed(42)
fa.cf <- cfQ(loadings(fa.unrotated), kappa = 0.3, normalize = TRUE,
             randomStarts = 100)
fa.cf

if (getRversion() >= "4.5.1") {
  # Single-step via factanal (correct in R >= 4.5.1)
  set.seed(42)
  fa.factanal <- factanal(factors = 3, covmat = NetherlandsTV,
                          normalize = TRUE, rotation = "cfQ",
                          control = list(rotate = list(kappa = 0.3,
                                                       randomStarts = 100)))
  fa.sorted <- print(fa.cf, sortLoadings = TRUE)
  cat("Maximum difference in loadings:\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")
}


###################################################
### code chunk number 26: GPA1guide.Stex:874-952
###################################################
  factanal_fit <- function(fa, R_obs, n) {
  # Compute common factor model fit indices from factanal output.
  # For MLE extraction only (factanal). For other estimators (minres,
  # principal axis, WLS) the chi-square statistic is not defined in
  # the same way and this function should not be used.
  #
  # Args:
  #   fa    : a factanal object (rotation = "none" recommended)
  #   R_obs : observed correlation or covariance matrix passed to factanal
  #   n     : sample size (n.obs passed to factanal)

  if (is.null(fa$STATISTIC))
    stop("factanal did not compute a chi-square statistic. ",
         "Ensure n.obs is specified when passing a covariance matrix.")

  # Ensure we work with a correlation matrix
  R_obs <- cov2cor(as.matrix(R_obs))
  p     <- nrow(R_obs) 
  L     <- loadings(fa)
  k     <- ncol(L)           # number of factors extracted
 
  # Model-implied correlation matrix using factanal uniquenesses
  # fa$uniquenesses are the MLE estimated unique variances
  R_hat <- L %*% t(L) + diag(fa$uniquenesses)
  R_hat <- cov2cor(R_hat)
  
  # Residual matrix --- off-diagonal only
  Resid <- R_obs - R_hat
  diag(Resid) <- 0

  # Test of the hypothesis that k factors are sufficient. Identical to factanal
  F_ml <- log(det(R_hat)) - log(det(R_obs)) + sum(diag(R_obs %*% solve(R_hat))) - p
  chi2 <- (n - 1 - (2 * p + 5)/6 - (2 * k)/3) * F_ml  
  df <- ((p - k)^2 - (p + k)) / 2
  pval <- pchisq(chi2, df, lower.tail = FALSE)

  # SRMR: standardized root mean square residual
  rstar.off <- sum(Resid^2) / 2
  srmr <- sqrt(rstar.off / (p * (p - 1)))

  # RMSEA: root mean square error of approximation
  rmsea <- sqrt(max(0, chi2 / (df * n) - 1 / (n - 1)))

  # Null model: all items uncorrelated
  chi2_null <- (n - 1) * sum(R_obs[lower.tri(R_obs)]^2)
  df_null   <- p * (p - 1) / 2

  # CFI: comparative fit index
  cfi <- (max(chi2_null - df_null, 0) - max(chi2 - df, 0)) /
          max(chi2_null - df_null, 0)

  # TLI: Tucker-Lewis index
  tli <- (chi2_null / df_null - chi2 / df) /
         (chi2_null / df_null - 1)

  # AIC and BIC
  aic <- chi2 - 2 * df
  bic <- chi2 - log(n) * df

  # Print results
  cat("Factor Model Fit Indices (MLE only)\n")
  cat("------------------------------------\n")
  cat(sprintf("Chi-square (df = %d):  %.3f  p = %.4f\n", df, chi2, pval))
  cat(sprintf("RMSEA:                 %.4f\n", rmsea))
  cat(sprintf("SRMR:                  %.4f\n", srmr))
  cat(sprintf("CFI:                   %.4f\n", cfi))
  cat(sprintf("TLI:                   %.4f\n", tli))
  cat(sprintf("AIC:                   %.3f\n", aic))
  cat(sprintf("BIC:                   %.3f\n", bic))
  cat("\nTop 5 absolute residuals:\n")
  resid_vals <- Resid[lower.tri(Resid)]
  print(round(sort(abs(resid_vals), decreasing = TRUE)[1:5], 4))

  invisible(c(chi2  = chi2,  df   = df,   pval = pval,
              rmsea = rmsea, srmr = srmr,
              cfi   = cfi,   tli  = tli,
              aic   = aic,   bic  = bic))
}


###################################################
### code chunk number 27: GPA1guide.Stex:958-966
###################################################
data("WansbeekMeijer", package = "GPArotation")
fa2 <- factanal(factors = 2, covmat = NetherlandsTV, rotation = "none")
fa3 <- factanal(factors = 3, covmat = NetherlandsTV, rotation = "none")
cat("=== 2 factors ===\n")
fit2 <- factanal_fit(fa2, cov2cor(NetherlandsTV$cov), n = 2154)

cat("\n=== 3 factors ===\n")
fit3 <- factanal_fit(fa3, cov2cor(NetherlandsTV$cov), n = 2154)


###################################################
### code chunk number 28: GPA1guide.Stex:969-977
###################################################
cat("Model comparison:\n")
cat(sprintf("           2-factor  3-factor\n"))
cat(sprintf("RMSEA:     %8.4f  %8.4f\n", fit2["rmsea"],  fit3["rmsea"]))
cat(sprintf("SRMR:      %8.4f  %8.4f\n", fit2["srmr"],   fit3["srmr"]))
cat(sprintf("CFI:       %8.4f  %8.4f\n", fit2["cfi"],    fit3["cfi"]))
cat(sprintf("TLI:       %8.4f  %8.4f\n", fit2["tli"],    fit3["tli"]))
cat(sprintf("AIC:       %8.2f  %8.2f\n", fit2["aic"],    fit3["aic"]))
cat(sprintf("BIC:       %8.2f  %8.2f\n", fit2["bic"],    fit3["bic"]))