fix factanal() internal sortLoading() needs to adapt rot.matrix, too

git-svn-id: https://svn.r-project.org/R/trunk@88047 00db46b3-68df-0310-9c12-caf00c1e9a41

Author: maechler

doc/NEWS.Rd

@@ -47,7 +47,7 @@
       results in confusing-to-read code.
     }
   }
-  
+
   \subsection{BUG FIXES}{
     \itemize{
       \item Setting \code{\link{attributes}} on primitive functions is an
@@ -810,6 +810,10 @@
       \item \command{R CMD check} now also reports bad symbols in
       package shared objects linked in from local static libraries
       (\PR{18789}).
+
+      \item \code{factanal()} now works correctly also, e.g., for
+      \pkg{GPArotation}, \code{oblimin()} rotations, fixing \PR{18417},
+      thanks to \I{Coen Bernaards} and others.
     }
   }
 }

src/library/stats/R/factanal.R

@@ -1,7 +1,7 @@
 #  File src/library/stats/R/factanal.R
 #  Part of the R package, https://www.R-project.org
 #
-#  Copyright (C) 1995-2015 The R Core Team
+#  Copyright (C) 1995-2025 The R Core Team
 #
 #  This program is free software; you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
@@ -16,31 +16,34 @@
 #  A copy of the GNU General Public License is available at
 #  https://www.R-project.org/Licenses/
 
-## Hmm, MM thinks diag(.) needs checking { diag(vec) when length(vec)==1 !}
-## However, MM does not understand that factor analysis
-##   is a *multi*variate technique!
 factanal <-
     function (x, factors, data = NULL, covmat = NULL, n.obs = NA,
               subset, na.action, start = NULL,
               scores = c("none", "regression", "Bartlett"),
               rotation = "varimax",
               control = NULL, ...)
 {
-    sortLoadings <- function(Lambda)
+    sortLoadings <- function(Lambda, rotm = NULL)
     {
         cn <- colnames(Lambda)
         Phi <- attr(Lambda, "covariance")
-        ssq <- apply(Lambda, 2L, function(x) -sum(x^2))
-        Lambda <- Lambda[, order(ssq), drop = FALSE]
-        colnames(Lambda) <- cn
+        io.ssq <- order(- if(is.null(Phi)) colSums(Lambda^2) else diag(Phi %*% crossprod(Lambda)))
+        Lambda <- Lambda[, io.ssq, drop = FALSE]
         neg <- colSums(Lambda) < 0
         Lambda[, neg] <- -Lambda[, neg]
-        if(!is.null(Phi)) {
-            unit <- ifelse(neg, -1, 1)
+        colnames(Lambda) <- cn
+        has.rot <- !is.null(rotm)
+        has.Ph  <- !is.null(Phi)
+        if(has.rot || has.Ph) # used for both
+            unit <- c(1, -1)[neg+1L] # = ifelse(neg, -1, 1)
+        ## FIXME:  <M>  %*%  diag(unit)  can be made faster (for non-small cases) below
+        if(has.rot)
+            rotm <- rotm %*% diag(unit)[, io.ssq]
+        if (has.Ph) {
             attr(Lambda, "covariance") <-
-                unit %*% Phi[order(ssq), order(ssq)] %*% unit
+                diag(unit) %*% Phi[io.ssq, io.ssq] %*% diag(unit)
         }
-        Lambda
+        list(Lambda = Lambda, rotm = rotm)
     }
     cl <- match.call()
     na.act <- NULL
@@ -98,6 +101,7 @@ factanal <-
                               "%d factor is too many for %d variables",
                               "%d factors are too many for %d variables"),
                      factors, p), domain = NA)
+    ## FIXME: cov2cor() shows how to do this faster (for large dim)
     sds <- sqrt(diag(cv))
     cv <- cv/(sds %o% sds)
 
@@ -121,8 +125,8 @@ factanal <-
     if(nc < 1) stop("no starting values supplied")
     best <- Inf
     for (i in 1L:nc) {
-        nfit <- factanal.fit.mle(cv, factors, start[, i],
-                                 max(cn$lower, 0), cn$opt)
+        ## factanal.fit.mle(cmat, factors, start=NULL, lower = 0.005, control = NULL, ...)
+        nfit <- factanal.fit.mle(cv, factors, start[, i], max(cn$lower, 0), cn$opt)
         if(cn$trace)
             cat("start", i, "value:", format(nfit$criteria[1L]),
                 "uniqs:", format(as.vector(round(nfit$uniquenesses, 4))), "\n")
@@ -136,20 +140,25 @@ factanal <-
                       "unable to optimize from this starting value",
                       "unable to optimize from these starting values"),
              domain = NA)
-    load <- fit$loadings
-    if(rotation != "none") {
+                                        # the following line changed by C. Bernaards, 20 December 2022
+    load <- sortLoadings(fit$loadings)$Lambda
+    if(is.function(rotation) || rotation != "none") {
         rot <- do.call(rotation, c(list(load), cn$rotate))
         load <- if (is.list(rot)) {
-          load <- rot$loadings
-          fit$rotmat <-
-              if(inherits(rot, "GPArotation")) t(solve(rot$Th))
-              else rot$rotmat
-          rot$loadings
-        } else rot
+                    fit$rotmat <-
+                        if(inherits(rot, "GPArotation")) t(solve(rot$Th))
+                        else rot$rotmat
+                    rot$loadings
+                } else rot
+        if(is.list(rot) && !is.null(rot$Phi)) attr(load, "covariance") <- rot$Phi
     }
-    fit$loadings <- sortLoadings(load)
+                                        # the following lines added by C. Bernaards, 20 December 2022
+    loadrot <- sortLoadings(load, fit$rotmat)
+    fit$loadings <- loadrot$Lambda
+    if(!is.null(loadrot$rotm)) fit$rotmat <- loadrot$rotm
+                                        # end additions C. Bernaards, 20 December 2022
     class(fit$loadings) <- "loadings"
-    fit$na.action <- na.act # not used currently
+    fit$na.action <- na.act # book keeping
     if(have.x && scores != "none") {
         Lambda <- fit$loadings
         zz <- scale(z, TRUE, TRUE)
@@ -217,8 +226,8 @@ factanal.fit.mle <-
         start <- (1 - 0.5*factors/p)/diag(solve(cmat))
     res <- optim(start, FAfn, FAgr, method = "L-BFGS-B",
                  lower = lower, upper = 1,
-                 control = c(list(fnscale=1,
-                 parscale = rep(0.01, length(start))), control),
+                 control = c(list(fnscale = 1,
+                                  parscale = rep(0.01, length(start))), control),
                  q = factors, S = cmat)
     Lambda <- FAout(res$par, cmat, factors)
     dimnames(Lambda) <- list(dimnames(cmat)[[1L]],
@@ -248,11 +257,13 @@ print.loadings <- function(x, digits = 3L, cutoff = 0.1, sort = FALSE, ...)
         Lambda <- Lambda[order(mx, 1L:p), , drop=FALSE]
     }
     cat("\nLoadings:\n")
-    fx <- setNames(format(round(Lambda, digits)), NULL)
+    fx <- setNames(format(round(Lambda, digits)), NULL) # char matrix
     nc <- nchar(fx[1L], type="c")
     fx[abs(Lambda) < cutoff] <- strrep(" ", nc)
     print(fx, quote = FALSE, ...)
-    vx <- colSums(x^2)
+    Phi <- attr(Lambda, "covariance")
+    vx <- if(is.null(Phi)) colSums(Lambda^2) else diag(Phi %*% crossprod(Lambda))
+    names(vx) <- colnames(Lambda)
     varex <- rbind("SS loadings" = vx)
     if(is.null(attr(x, "covariance"))) {
         varex <- rbind(varex, "Proportion Var" = vx/p)
@@ -269,19 +280,19 @@ print.factanal <- function(x, digits = 3, ...)
     cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
     cat("Uniquenesses:\n")
     print(round(x$uniquenesses, digits), ...)
-    print(x$loadings, digits = digits, ...)
+    print(x$loadings, digits = digits, ...) #-> print.loadings()
                                         # the following lines added by J. Fox, 26 June 2005
-    if (!is.null(x$rotmat)){
-      tmat <- solve(x$rotmat)
-      R <- tmat %*% t(tmat)
-      factors <- x$factors
-      rownames(R) <- colnames(R) <- paste0("Factor", 1:factors)
+    if(!is.null(x$rotmat)) {
+        tmat <- solve(x$rotmat)
+        R <- tcrossprod(tmat) # t . t'
+        factors <- x$factors
+        rownames(R) <- colnames(R) <- paste0("Factor", 1:factors)
 
                                         # the following line changed by Ulrich Keller, 9 Sept 2008
-      if (!isTRUE(all.equal(c(R), c(diag(factors))))) {
-        cat("\nFactor Correlations:\n")
-        print(R, digits=digits, ...)
-      }
+        if (!isTRUE(all.equal(c(R), c(diag(factors))))) {
+            cat("\nFactor Correlations:\n")
+            print(R, digits=digits, ...)
+        }
     }
                                         # end additions J. Fox, 23 June 2005
     if(!is.null(x$STATISTIC)) {

src/library/stats/man/factanal.Rd

@@ -1,6 +1,6 @@
 % File src/library/stats/man/factanal.Rd
 % Part of the R package, https://www.R-project.org
-% Copyright 1995-2018 R Core Team
+% Copyright 1995-2025 R Core Team
 % Distributed under GPL 2 or later
 
 \name{factanal}
@@ -18,44 +18,46 @@ factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
          rotation = "varimax", control = NULL, \dots)
 }
 \arguments{
-  \item{x}{A formula or a numeric matrix or an object that can be
+  \item{x}{\code{\link{formula}} or a numeric matrix or an object that can be
     coerced to a numeric matrix.}
-  \item{factors}{The number of factors to be fitted.}
-  \item{data}{An optional data frame (or similar: see
+  \item{factors}{integer number of factors to be fitted.}
+  \item{data}{an optional data frame (or similar: see
     \code{\link{model.frame}}), used only if \code{x} is a formula.  By
     default the variables are taken from \code{environment(formula)}.}
-  \item{covmat}{A covariance matrix, or a covariance list as returned by
+  \item{covmat}{covariance matrix or a covariance list as returned by
     \code{\link{cov.wt}}.  Of course, correlation matrices are covariance
     matrices.}
-  \item{n.obs}{The number of observations, used if \code{covmat} is a
+  \item{n.obs}{the number of observations, used if \code{covmat} is a
     covariance matrix.}
-  \item{subset}{A specification of the cases to be used, if \code{x} is
+  \item{subset}{a specification of the cases to be used, if \code{x} is
     used as a matrix or formula.}
-  \item{na.action}{The \code{na.action} to be used if \code{x} is
-    used as a formula.}
+  \item{na.action}{the \code{\link{na.action}} to be used if \code{x} is
+    used as a formula, in which case \code{\link{model.frame}()} is
+    called.}
   \item{start}{\code{NULL} or a matrix of starting values, each column
     giving an initial set of uniquenesses.}
   \item{scores}{Type of scores to produce, if any.  The default is none,
     \code{"regression"} gives Thompson's scores, \code{"Bartlett"} given
     Bartlett's weighted least-squares scores. Partial matching allows
     these names to be abbreviated.}
-  \item{rotation}{character. \code{"none"} or the name of a function
+  \item{rotation}{a \code{\link{function}} or \code{\link{character}}
+    string.  In the latter case, \code{"none"} or the name of a function
     to be used to rotate the factors: it will be called with first
     argument the loadings matrix, and should return a list with component
     \code{loadings} giving the rotated loadings, or just the rotated loadings.}
-  \item{control}{A list of control values,
+  \item{control}{a \code{\link{list}} of control values with components
     \describe{
-      \item{nstart}{The number of starting values to be tried if
+      \item{\code{nstart}}{The number of starting values to be tried if
         \code{start = NULL}. Default 1.}
-      \item{trace}{logical. Output tracing information? Default \code{FALSE}.}
-      \item{lower}{The lower bound for uniquenesses during
+      \item{\code{trace}}{logical. Output tracing information? Default \code{FALSE}.}
+      \item{\code{lower}}{The lower bound for uniquenesses during
         optimization. Should be > 0. Default 0.005.}
-      \item{opt}{A list of control values to be passed to
-        \code{\link{optim}}'s \code{control} argument.}
-      \item{rotate}{a list of additional arguments for the rotation function.}
+      \item{\code{opt}}{a list of control values to be passed to
+	\code{\link{optim}}'s \code{control} argument.}
+      \item{\code{rotate}}{a list of additional arguments for the rotation function.}
     }
   }
-  \item{\dots}{Components of \code{control} can also be supplied as
+  \item{\dots}{components of \code{control} can also be supplied directly as
     named arguments to \code{factanal}.}
 }
 \details{
@@ -128,7 +130,8 @@ factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
   \item{loadings}{A matrix of loadings, one column for each factor.  The
     factors are ordered in decreasing order of sums of squares of
     loadings, and given the sign that will make the sum of the loadings
-    positive.  This is of class \code{"loadings"}: see
+    positive. For correlated factors, the correlation matrix is
+    stored as an attribute labeled "covariance". This is of class \code{"loadings"}: see
     \code{\link{loadings}} for its \code{print} method.}
   \item{uniquenesses}{The uniquenesses computed.}
   \item{correlation}{The correlation matrix used.}

src/library/stats/man/loadings.Rd

@@ -1,19 +1,20 @@
 % File src/library/stats/man/loadings.Rd
 % Part of the R package, https://www.R-project.org
-% Copyright 1995-2014 R Core Team
+% Copyright 1995-2025 R Core Team
 % Distributed under GPL 2 or later
 
 \name{loadings}
+\title{Print Factor Analysis 'factanal()', Notably Loadings}
 \alias{loadings}
 \alias{print.loadings}
 \alias{print.factanal}
-\title{Print Loadings in Factor Analysis}
 \description{
-  Extract or print loadings in factor analysis (or principal
-  components analysis).
+  Print factor analysis result,  \code{\link{factanal}()};
+  extract or print loadings in factor analysis (or principal
+  components analysis, PCA).
 }
 \usage{
-loadings(x, ...)
+loadings(x, \dots)
 
 \method{print}{loadings}(x, digits = 3, cutoff = 0.1, sort = FALSE, \dots)
 
@@ -31,8 +32,9 @@ loadings(x, ...)
     than 0.5 (in modulus) is assigned to the factor with the largest
     loading, and the variables are printed in the order of the factor
     they are assigned to, then those unassigned.}
-  \item{\dots}{further arguments for other methods,
-    ignored for \code{loadings}.}
+  \item{\dots}{further arguments for other methods; for \code{factanal}
+    printing, notably passed to \code{print(<loadings>, ...)}, hence, e.g.,
+    \code{cutoff = <frac>}; ignored for \code{loadings}.}
 }
 
 \details{
@@ -47,7 +49,7 @@ loadings(x, ...)
   The \code{print} method for class \code{"\link{factanal}"} calls the
   \code{"loadings"} method to print the loadings, and so passes down
   arguments such as \code{cutoff} and \code{sort}.
-  
+
   The signs of the loadings vectors are arbitrary for both factor
   analysis and PCA.
 }

src/library/stats/tests/factanal-tst.R

@@ -0,0 +1,42 @@
+## Subject: [Bug 18417] factanal returns useless SS loading values
+## Date: Sat, 22 Mar 2025 23:21:10 +0000
+## --- Comment #21 from Coen Bernaards
+
+## create unrotated loadings matrix
+data(ability.cov, package="datasets")
+res <- factanal(covmat = ability.cov$cov, factors = 3, rotation = "none")
+
+if(FALSE) {
+    ## create the output from oblimin of GPArotation as a list;
+    require(GPArotation)
+    set.seed(100)
+    obli100 <- oblimin(res$loadings, Tmat = Random.Start(3), eps = 1e-8, maxit = 10000)
+    ## keep only what's needed:
+    obli100[setdiff(names(obli100),c("loadings", "Phi", "Th", "orthogonal"))] <- NULL
+    saveRDS(obli100, "obli100.rds")
+}
+
+obli100 <- readRDS("obli100.rds") # in ../stats/tests/
+if(getRversion() >= "4.6") { ## perform the factanal using this oblimin rotation [ rotation now may be _function_
+  resO100 <- factanal(covmat = ability.cov$cov, factors = 3, rotation = \(...) obli100)
+} else { ## old R:
+  obli_FN <- function(...) obli100
+  resO100 <- factanal(covmat = ability.cov$cov, factors = 3, rotation = "obli_FN")
+}
+resO100 # print.factanal(..)  --> print.loadings(..)
+
+oriL <-   res     $ loadings
+ldns   <- resO100 $ loadings
+rotmat <- resO100 $ rotmat
+# no longer can the unrotated matrix be recreated. the max abs diff between
+# unrotated and recreated unrotated loadings should be near zero but is not:
+print.default(dev <- ldns %*% solve(rotmat) - oriL)
+stopifnot(max(abs(dev)) < 1e-14) # max|dev| was 2.000982 (R <= 4.4.x), now 2.22e-16 [Lnx x86_64, gcc]
+
+ldns # prints loadings too
+## check the "SS loadings .... " line
+L <- ldns
+Phi <- attr(L, "covariance")
+(SSld <- diag(Phi %*% crossprod(L)))
+##      1.896492 1.124306 1.037336
+stopifnot(all.equal(c(1.8964924541, 1.1243062406, 1.0373362639), SSld)) # seen  8.49e-12