Nothing
# various rotation criteria and their gradients
# YR 05 April 2019: initial version
# YR 14 June 2019: add more rotation criteria
# references:
#
# Bernaards, C. A., & Jennrich, R. I. (2005). Gradient projection algorithms
# and software for arbitrary rotation criteria in factor analysis. Educational
# and Psychological Measurement, 65(5), 676-696.
# old website: http://web.archive.org/web/20180708170331/http://www.stat.ucla.edu/research/gpa/splusfunctions.net
#
# Browne, M. W. (2001). An overview of analytic rotation in exploratory factor
# analysis. Multivariate behavioral research, 36(1), 111-150.
#
# Mulaik, S. A. (2010). Foundations of factor analysis (Second Edition).
# Boca Raton: Chapman and Hall/CRC.
# Note: this is YR's implementation, not a copy of the GPArotation
# package
#
# Why did I write my own functions (and not use the GPArotation):
# - to better understand what is going on
# - to have direct access to the gradient functions
# - to avoid yet another dependency
# - to simplify further experiments
# Orthomax family (Harman, 1960)
#
# gamma = 0 -> quartimax
# gamma = 1/2 -> biquartimax
# gamma = 1/P -> equamax
# gamma = 1 -> varimax
#
lav_matrix_rotate_orthomax <- function(LAMBDA = NULL, orthomax.gamma = 1,
..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
# center L2 column-wise
cL2 <- t(t(L2) - orthomax.gamma * colMeans(L2))
out <- -1 * sum(L2 * cL2) / 4
if (grad) {
attr(out, "grad") <- -1 * LAMBDA * cL2
}
out
}
# Crawford-Ferguson (1970) family
#
# combine penalization for 1) row complexity, and 2) column complexity
# if combined with orthogonal rotation, this is equivalent to the
# orthomax family:
#
# quartimax -> gamma = 0 (only row complexity)
# varimax -> gamma = 1/nrow
# equamax -> gamma = ncol/(2*nrow)
# parsimax -> gamma = (ncol - 1)/(nrow + ncol - 2)
# factor parsimony -> gamma = 1 (only column complexity)
#
# the Crawford-Ferguson family is also equivalent to the oblimin family
# if the latter is restricted to orthogonal rotation
#
lav_matrix_rotate_cf <- function(LAMBDA = NULL, cf.gamma = 0, ...,
grad = FALSE) {
# check if gamma is between 0 and 1?
nRow <- nrow(LAMBDA)
nCol <- ncol(LAMBDA)
ROW1 <- matrix(1.0, nCol, nCol)
diag(ROW1) <- 0.0
COL1 <- matrix(1.0, nRow, nRow)
diag(COL1) <- 0.0
L2 <- LAMBDA * LAMBDA
LR <- L2 %*% ROW1
LC <- COL1 %*% L2
f1 <- sum(L2 * LR) / 4
f2 <- sum(L2 * LC) / 4
out <- (1 - cf.gamma) * f1 + cf.gamma * f2
if (grad) {
attr(out, "grad") <- ((1 - cf.gamma) * LAMBDA * LR) +
(cf.gamma * LAMBDA * LC)
}
out
}
# Oblimin family (Carroll, 1960; Harman, 1976)
#
# quartimin -> gamma = 0
# biquartimin -> gamma = 1/2
# covarimin -> gamma = 1
#
# if combined with orthogonal rotation, this is equivalent to the
# orthomax family (they have the same optimizers):
#
# gamma = 0 -> quartimax
# gamma = 1/2 -> biquartimax
# gamma = 1 -> varimax
# gamma = P/2 -> equamax
#
lav_matrix_rotate_oblimin <- function(LAMBDA = NULL, oblimin.gamma = 0, ...,
grad = FALSE) {
nRow <- nrow(LAMBDA)
nCol <- ncol(LAMBDA)
ROW1 <- matrix(1.0, nCol, nCol)
diag(ROW1) <- 0.0
L2 <- LAMBDA * LAMBDA
LR <- L2 %*% ROW1
Jp <- matrix(1, nRow, nRow) / nRow
# see Jennrich (2002, p. 11)
tmp <- (diag(nRow) - oblimin.gamma * Jp) %*% LR
# same as t( t(L2) - gamma * colMeans(L2) ) %*% ROW1
out <- sum(L2 * tmp) / 4
if (grad) {
attr(out, "grad") <- LAMBDA * tmp
}
out
}
# quartimax criterion
# Carroll (1953); Saunders (1953) Neuhaus & Wrigley (1954); Ferguson (1954)
# we use here the equivalent 'Ferguson, 1954' variant
# (See Mulaik 2010, p. 303)
lav_matrix_rotate_quartimax <- function(LAMBDA = NULL, ..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
out <- -1 * sum(L2 * L2) / 4
if (grad) {
attr(out, "grad") <- -1 * LAMBDA * L2
}
out
}
# varimax criterion
# Kaiser (1958, 1959)
#
# special case of the Orthomax family (Harman, 1960), where gamma = 1
# see Jennrich (2001, p. 296)
lav_matrix_rotate_varimax <- function(LAMBDA = NULL, ..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
# center L2 column-wise
cL2 <- t(t(L2) - colMeans(L2))
out <- -1 * abs(sum(L2 * cL2)) / 4 # abs needed?
if (grad) {
attr(out, "grad") <- -1 * LAMBDA * cL2
}
out
}
# quartimin criterion (part of Carroll's oblimin family
lav_matrix_rotate_quartimin <- function(LAMBDA = NULL, ..., grad = FALSE) {
nCol <- ncol(LAMBDA)
ROW1 <- matrix(1.0, nCol, nCol)
diag(ROW1) <- 0.0
L2 <- LAMBDA * LAMBDA
LR <- L2 %*% ROW1
out <- sum(L2 * LR) / 4
if (grad) {
attr(out, "grad") <- LAMBDA * LR
}
out
}
# Browne's (2001) version of Yates (1984) geomin criterion
#
# we use the exp/log trick as in Bernaard & Jennrich (2005, p. 687)
lav_matrix_rotate_geomin <- function(LAMBDA = NULL, geomin.epsilon = 0.01,
..., grad = FALSE) {
nCol <- ncol(LAMBDA)
L2 <- LAMBDA * LAMBDA
L2 <- L2 + geomin.epsilon
if (geomin.epsilon < sqrt(.Machine$double.eps)) {
# Yates's original formula
tmp <- apply(L2, 1, prod)^(1 / nCol)
} else {
tmp <- exp(rowSums(log(L2)) / nCol)
}
out <- sum(tmp)
if (grad) {
attr(out, "grad") <- (2 / nCol) * LAMBDA / L2 * tmp
}
out
}
# simple entropy
# seems to only work for orthogonal rotation
lav_matrix_rotate_entropy <- function(LAMBDA = NULL, ..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
# handle zero elements -> replace by '1', so log(1) == 0
L2[L2 == 0] <- 1
out <- -1 * sum(L2 * log(L2)) / 2
if (grad) {
attr(out, "grad") <- -LAMBDA * log(L2) - LAMBDA
}
out
}
# McCammon's (1966) Minimum Entropy Criterion
#
# for p-vector x, where x > 0 and sum(x) = 1, we have
# - entropy(x) == 0, if there is only one 1, and all zeroes
# - entropy(x) == max == log(p) if all elements are 1/p
# - entropy(x) is similar as complexity(x), but also measure of equality
# of elements of x
#
# works only ok with orthogonal rotation!
lav_matrix_rotate_mccammon <- function(LAMBDA = NULL, ..., grad = FALSE) {
nCol <- ncol(LAMBDA)
nRow <- nrow(LAMBDA)
L2 <- LAMBDA * LAMBDA
# entropy function (Browne, 2001, eq 9)
f_entropy <- function(x) {
-1 * sum(ifelse(x > 0, x * log(x), 0))
}
# sums of rows/columns/all
sumi. <- rowSums(L2)
sum.j <- colSums(L2)
sum.. <- sum(L2)
Q1 <- f_entropy(t(L2) / sum.j) # encouraging columns with few large,
# and many small elements
Q2 <- f_entropy(sum.j / sum..) # encouraging equal column sums
# minimize
out <- log(Q1) - log(Q2)
if (grad) { # See Bernaards and Jennrich 2005 page 685+686
H <- -(log(t(t(L2) / sum.j)) + 1)
G1 <- t(t(H) / sum.j - rowSums(t(L2 * H) / (sum.j * sum.j)))
h <- -(log(sum.j / sum..) + 1)
alpha <- as.numeric(h %*% sum.j) / (sum.. * sum..) # paper divides by
# sum.., not sum..^2??
G2 <- matrix(h / sum.. - alpha, nRow, nCol, byrow = TRUE)
attr(out, "grad") <- 2 * LAMBDA * (G1 / Q1 - G2 / Q2)
}
out
}
# Infomax
# McKeon (1968, unpublished) and Browne (2001)
# Treat LAMBDA^2 as a contingency table, and use simplicity function based
# on tests for association; most effective was LRT for association
# (see Agresti, 1990, eq 3.13) which is maximized for max simplicity
#
# McKeon: criterion may be regarded as a measure of information about row
# categories conveyed by column categories (and vice versa); hence infomax
# - favors perfect cluster
# - discourages general factor
# - both for orthogonal and oblique rotation
#
# Note: typo in Browne (2001), see last paragraph of Bernaards and
# Jennrich (2005) page 684
lav_matrix_rotate_infomax <- function(LAMBDA = NULL, ..., grad = FALSE) {
nCol <- ncol(LAMBDA)
nRow <- nrow(LAMBDA)
L2 <- LAMBDA * LAMBDA
# entropy function (Browne, 2001, eq 9)
f_entropy <- function(x) {
-1 * sum(ifelse(x > 0, x * log(x), 0))
}
# sums of rows/columns/all
sumi. <- rowSums(L2)
sum.j <- colSums(L2)
sum.. <- sum(L2)
Q1 <- f_entropy(L2 / sum..) # Bernaards & Jennrich version!! (Browne
# divides by sum.j, like in McCammon)
Q2 <- f_entropy(sum.j / sum..)
Q3 <- f_entropy(sumi. / sum..)
# minimize
out <- log(nCol) + Q1 - Q2 - Q3
if (grad) {
H <- -(log(L2 / sum..) + 1)
alpha <- sum(L2 * H) / (sum.. * sum..)
G1 <- H / sum.. - alpha
hj <- -(log(sum.j / sum..) + 1)
alphaj <- as.numeric(hj %*% sum.j) / (sum.. * sum..)
G2 <- matrix(hj, nRow, nCol, byrow = TRUE) / sum.. - alphaj
hi <- -(log(sumi. / sum..) + 1)
alphai <- as.numeric(sumi. %*% hi) / (sum.. * sum..)
G3 <- matrix(hi, nRow, nCol) / sum.. - alphai
attr(out, "grad") <- 2 * LAMBDA * (G1 - G2 - G3)
}
out
}
# oblimax
# Harman, 1976; Saunders, 1961
#
# for orthogonal rotation, oblimax is equivalent to quartimax
lav_matrix_rotate_oblimax <- function(LAMBDA = NULL, ..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
# minimize version
out <- -log(sum(L2 * L2)) + 2 * log(sum(L2))
if (grad) {
attr(out, "grad") <- (-4 * L2 * LAMBDA / (sum(L2 * L2))
+ 4 * LAMBDA / (sum(L2)))
}
out
}
# Bentler's Invariant Pattern Simplicity
# Bentler (1977)
#
#
lav_matrix_rotate_bentler <- function(LAMBDA = NULL, ..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
L2tL2 <- crossprod(L2)
L2tL2.inv <- lav_matrix_symmetric_inverse(S = L2tL2, logdet = TRUE)
L2tL2.logdet <- attr(L2tL2.inv, "logdet")
DIag <- diag(L2tL2)
DIag.inv <- diag(1 / DIag)
DIag.logdet <- sum(log(DIag)) # add small constant?
# minimize version
out <- -(L2tL2.logdet - DIag.logdet) / 4
if (grad) {
attr(out, "grad") <- -LAMBDA * (L2 %*% (L2tL2.inv - DIag.inv))
}
out
}
# The Tandem criteria
# Comrey (1967)
#
# only for sequential use:
# - tandem1 is used to determine the number of factors
# (it removes the minor factors)
# - tandomII is used for final rotation
#
lav_matrix_rotate_tandem1 <- function(LAMBDA, ..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
LL <- tcrossprod(LAMBDA)
LL2 <- LL * LL
# minimize version
out <- -1 * sum(L2 * (LL2 %*% L2))
if (grad) {
tmp1 <- 4 * LAMBDA * (LL2 %*% L2)
tmp2 <- 4 * (LL * (L2 %*% t(L2))) %*% LAMBDA
attr(out, "grad") <- -tmp1 - tmp2
}
out
}
lav_matrix_rotate_tandem2 <- function(LAMBDA, ..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
LL <- tcrossprod(LAMBDA)
LL2 <- LL * LL
# minimize version
out <- sum(L2 * ((1 - LL2) %*% L2))
if (grad) {
tmp1 <- 4 * LAMBDA * ((1 - LL2) %*% L2)
tmp2 <- 4 * (LL * tcrossprod(L2, L2)) %*% LAMBDA
attr(out, "grad") <- tmp1 - tmp2
}
out
}
# simplimax
# Kiers (1994)
#
# oblique rotation method
# designed to rotate so that a given number 'k' of small loadings are
# as close to zero as possible
#
# may be viewed as partially specified target rotation with
# dynamically chosen weights
#
lav_matrix_rotate_simplimax <- function(LAMBDA = NULL, k = nrow(LAMBDA),
..., grad = FALSE) {
L2 <- LAMBDA * LAMBDA
# 'k' smallest element of L2
small.element <- sort(L2)[k]
# which elements are smaller than (or equal than) 'small.element'?
ID <- sign(L2 <= small.element)
# minimize version
out <- sum(L2 * ID)
if (grad) {
attr(out, "grad") <- 2 * ID * LAMBDA
}
out
}
# target rotation
# Harman, 1976
#
# LAMBDA is rotated toward a specified target matrix 'target'
#
# Note: 'target' must be fully specified; if there are any NAs
# use lav_matrix_rotate_pst() instead
#
lav_matrix_rotate_target <- function(LAMBDA = NULL, target = NULL,
..., grad = FALSE) {
# squared difference
DIFF <- LAMBDA - target
DIFF2 <- DIFF * DIFF
out <- sum(DIFF2, na.rm = TRUE)
if (grad) {
tmp <- 2 * DIFF
# change NAs to zero
tmp[is.na(tmp)] <- 0
attr(out, "grad") <- tmp
}
out
}
# partially specified target rotation
#
# Browne 1972a, 1972b
#
# a pre-specified weight matrix W with ones/zeroes determines
# which elements of (LAMBDA - target) are used by the rotation criterion
#
# if 'target' contains NAs, they should correspond to '0' values in the
# target.mask matrix
#
lav_matrix_rotate_pst <- function(LAMBDA = NULL, target = NULL,
target.mask = NULL, ..., grad = FALSE) {
# mask target+LAMBDA
target <- target.mask * target
LAMBDA <- target.mask * LAMBDA
# squared difference
DIFF <- LAMBDA - target
DIFF2 <- DIFF * DIFF
# minimize
out <- sum(DIFF2, na.rm = TRUE)
if (grad) {
tmp <- 2 * DIFF
# change NAs to zero
tmp[is.na(tmp)] <- 0
attr(out, "grad") <- tmp
}
out
}
# bi-quartimin
#
# Jennrich & Bentler 2011
#
lav_matrix_rotate_biquartimin <- function(LAMBDA, ..., grad = FALSE) {
# see Matlab code page 549
stopifnot(ncol(LAMBDA) > 1L)
# remove first column
LAMBDA.group <- LAMBDA[, -1, drop = FALSE]
# apply quartimin on the 'group' part
out <- lav_matrix_rotate_quartimin(LAMBDA.group, ..., grad = grad)
if (grad) {
tmp <- attr(out, "grad")
attr(out, "grad") <- cbind(0, tmp)
}
out
}
# bi-geomin
#
# Jennrich & Bentler 2012
#
lav_matrix_rotate_bigeomin <- function(LAMBDA, geomin.epsilon = 0.01, ...,
grad = FALSE) {
stopifnot(ncol(LAMBDA) > 1L)
# remove first column
LAMBDA.group <- LAMBDA[, -1, drop = FALSE]
# apply geomin on the 'group' part
out <- lav_matrix_rotate_geomin(LAMBDA.group,
geomin.epsilon = geomin.epsilon, ...,
grad = grad
)
if (grad) {
tmp <- attr(out, "grad")
attr(out, "grad") <- cbind(0, tmp)
}
out
}
# gradient check
ilav_matrix_rotate_grad_test <- function(crit = NULL, ...,
LAMBDA = NULL,
nRow = 20L, nCol = 5L) {
# test matrix
if (is.null(LAMBDA)) {
LAMBDA <- matrix(rnorm(nRow * nCol), nRow, nCol)
}
ff <- function(x, ...) {
Lambda <- matrix(x, nRow, nCol)
crit(Lambda, ..., grad = FALSE)
}
GQ1 <- matrix(
numDeriv::grad(func = ff, x = as.vector(LAMBDA), ...),
nRow, nCol
)
GQ2 <- attr(crit(LAMBDA, ..., grad = TRUE), "grad")
if (lav_verbose()) {
print(list(LAMBDA = LAMBDA, GQ1 = GQ1, GQ2 = GQ2))
}
all.equal(GQ1, GQ2, tolerance = 1e-07)
}
ilav_matrix_rotate_grad_test_all <- function() {
# Orthomax family with various values for gamma
for (gamma in seq(0, 1, 0.2)) {
check <- ilav_matrix_rotate_grad_test(
crit = lav_matrix_rotate_orthomax,
gamma = gamma
)
if (is.logical(check) && check) {
cat(
"orthomax + gamma = ", sprintf("%3.1f", gamma),
": OK\n"
)
} else {
cat(
"orthomax + gamma = ", sprintf("%3.1f", gamma),
": FAILED\n"
)
}
}
# Crawford-Ferguson with various values for gamma
for (gamma in seq(0, 1, 0.2)) {
check <- ilav_matrix_rotate_grad_test(
crit = lav_matrix_rotate_cf,
gamma = gamma
)
if (is.logical(check) && check) {
cat(
"Crawford-Ferguson + gamma = ", sprintf("%3.1f", gamma),
": OK\n"
)
} else {
cat(
"Crawford-Ferguson + gamma = ", sprintf("%3.1f", gamma),
": FAILED\n"
)
}
}
# Oblimin family with various values for gamma
for (gamma in seq(0, 1, 0.2)) {
check <- ilav_matrix_rotate_grad_test(
crit = lav_matrix_rotate_oblimin,
gamma = gamma
)
if (is.logical(check) && check) {
cat(
"Oblimin + gamma = ", sprintf("%3.1f", gamma),
": OK\n"
)
} else {
cat(
"Oblimin + gamma = ", sprintf("%3.1f", gamma),
": FAILED\n"
)
}
}
# quartimax
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_quartimax)
if (is.logical(check) && check) {
cat("quartimax: OK\n")
} else {
cat("quartimax: FAILED\n")
}
# varimax
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_varimax)
if (is.logical(check) && check) {
cat("varimax: OK\n")
} else {
cat("varimax: FAILED\n")
}
# quartimin
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_quartimin)
if (is.logical(check) && check) {
cat("quartimin: OK\n")
} else {
cat("quartimin: FAILED\n")
}
# geomin
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_geomin)
if (is.logical(check) && check) {
cat("geomin: OK\n")
} else {
cat("geomin: FAILED\n")
}
# simple entropy
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_entropy)
if (is.logical(check) && check) {
cat("entropy: OK\n")
} else {
cat("entropy: FAILED\n")
}
# McCammon entropy criterion
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_mccammon)
if (is.logical(check) && check) {
cat("McCammon: OK\n")
} else {
cat("McCammon: FAILED\n")
}
# infomax
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_infomax)
if (is.logical(check) && check) {
cat("infomax: OK\n")
} else {
cat("infomax: FAILED\n")
}
# oblimax
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_oblimax)
if (is.logical(check) && check) {
cat("oblimax: OK\n")
} else {
cat("oblimax: FAILED\n")
}
# bentler
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_bentler)
if (is.logical(check) && check) {
cat("bentler: OK\n")
} else {
cat("bentler: FAILED\n")
}
# simplimax
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_simplimax)
if (is.logical(check) && check) {
cat("simplimax: OK\n")
} else {
cat("simplimax: FAILED\n")
}
# tandem1
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_tandem1)
if (is.logical(check) && check) {
cat("tandem1: OK\n")
} else {
cat("tandem1: FAILED\n")
}
# tandem2
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_tandem2)
if (is.logical(check) && check) {
cat("tandem2: OK\n")
} else {
cat("tandem2: FAILED\n")
}
# bi-quartimin
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_biquartimin)
if (is.logical(check) && check) {
cat("biquartimin: OK\n")
} else {
cat("biquartimin: FAILED\n")
}
# bi-quartimin
check <- ilav_matrix_rotate_grad_test(crit = lav_matrix_rotate_bigeomin)
if (is.logical(check) && check) {
cat("bigeomin: OK\n")
} else {
cat("bigeomin: FAILED\n")
}
}
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.