From jfox at mcmaster.ca  Tue Sep  9 16:53:30 2008
From: jfox at mcmaster.ca (John Fox)
Date: Tue, 9 Sep 2008 10:53:30 -0400
Subject: [Rd] Addendum to wishlist bug report #10931 (factanal)
	(PR#12754)
In-Reply-To: <20080909115511.4FF112834169@mail.pubhealth.ku.dk>
References: <20080909115511.4FF112834169@mail.pubhealth.ku.dk>
Message-ID: <000a01c9128b$d31dcef0$79596cd0$@ca>

Dear Ulrich,

I'd frankly forgotten about this, but can't see an argument for not making
this (or a similar) change. 

Thanks for reviving the issue.

John

------------------------------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox


> -----Original Message-----
> From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org]
On
> Behalf Of ulrich.keller at uni.lu
> Sent: September-09-08 7:55 AM
> To: r-devel at stat.math.ethz.ch
> Cc: R-bugs at r-project.org
> Subject: [Rd] Addendum to wishlist bug report #10931 (factanal) (PR#12754)
> 
> --=-hiYzUeWcRJ/+kx41aPIZ
> Content-Type: text/plain; charset="UTF-8"
> Content-Transfer-Encoding: 8bit
> 
> Hi,
> 
> on March 10 I filed a wishlist bug report asking for the inclusion of
> some changes to factanal() and the associated print method. The changes
> were originally proposed by John Fox in 2005; they make print.factanal()
> display factor correlations if factanal() is called with rotation =
> "promax". Since I got no replies, and I am really tired of my R-curious
> social science colleagues asking "What, it can't even display factor
> correlations?!", I made the changes myself. I would be very grateful if
> they'd find their way into a release.
> 
> I corrected a small error in John Fox's code and made another change
> that enables factor correlations not only for promax, but also for the
> rotation methods in package GPArotation.
> 
> The changes are against R-devel, downloaded on September 9th 2008.
> Changes are indicated by comments from John Fox and me. I also changed
> factanal.Rd accordingly, this is commented too.
> 
> My bug report is at
> http://bugs.r-project.org/cgi-bin/R/wishlist?id=10931;user=guest
> 
> John Fox's original post is at
> http://tolstoy.newcastle.edu.au/R/devel/05/06/1414.html
> 
> The changed files factanal.R and factanal.Rd are attached. If there is
> anything else I can do to help these changes make it into R, please let
> me know.
> 
> Thanks and best regards,
> 
> Uli
> 
> --
> Ulrich Keller
> Universit?? du Luxembourg
> EMACS research unit
> B.P. 2
> L-7201 Walferdange
> Luxembourg
> Mail ulrich.keller at uni.lu
> Phone +352 46 66 44 9 278
> 
> --=-hiYzUeWcRJ/+kx41aPIZ
> Content-Disposition: attachment; filename="factanal.R"
> Content-Type: text/plain; name="factanal.R"; charset="utf-8"
> Content-Transfer-Encoding: 7bit
> 
> #  File src/library/stats/R/factanal.R
> #  Part of the R package, http://www.R-project.org
> #
> #  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
> #  the Free Software Foundation; either version 2 of the License, or
> #  (at your option) any later version.
> #
> #  This program is distributed in the hope that it will be useful,
> #  but WITHOUT ANY WARRANTY; without even the implied warranty of
> #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
> #  GNU General Public License for more details.
> #
> #  A copy of the GNU General Public License is available at
> #  http://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)
>     {
>         cn <- colnames(Lambda)
>         Phi <- attr(Lambda, "covariance")
>         ssq <- apply(Lambda, 2, function(x) -sum(x^2))
>         Lambda <- Lambda[, order(ssq), drop = FALSE]
>         colnames(Lambda) <- cn
>         neg <- colSums(Lambda) < 0
>         Lambda[, neg] <- -Lambda[, neg]
>         if(!is.null(Phi)) {
>             unit <- ifelse(neg, -1, 1)
>             attr(Lambda, "covariance") <-
>                 unit %*% Phi[order(ssq), order(ssq)] %*% unit
>         }
>         Lambda
>     }
>     cl <- match.call()
>     na.act <- NULL
>     if (is.list(covmat)) {
>         if (any(is.na(match(c("cov", "n.obs"), names(covmat)))))
>             stop("'covmat' is not a valid covariance list")
>         cv <- covmat$cov
>         n.obs <- covmat$n.obs
>         have.x <- FALSE
>     }
>     else if (is.matrix(covmat)) {
>         cv <- covmat
>         have.x <- FALSE
>     }
>     else if (is.null(covmat)) {
>         if(missing(x)) stop("neither 'x' nor 'covmat' supplied")
>         have.x <- TRUE
>         if(inherits(x, "formula")) {
>             ## this is not a `standard' model-fitting function,
>             ## so no need to consider contrasts or levels
>             mt <- terms(x, data = data)
>             if(attr(mt, "response") > 0)
>                 stop("response not allowed in formula")
>             attr(mt, "intercept") <- 0
>             mf <- match.call(expand.dots = FALSE)
>             names(mf)[names(mf) == "x"] <- "formula"
>             mf$factors <- mf$covmat <- mf$scores <- mf$start <-
>                 mf$rotation <- mf$control <- mf$... <- NULL
>             mf[[1]] <- as.name("model.frame")
>             mf <- eval.parent(mf)
>             na.act <- attr(mf, "na.action")
>             if (.check_vars_numeric(mf))
>                 stop("factor analysis applies only to numerical
variables")
>             z <- model.matrix(mt, mf)
>         } else {
>             z <- as.matrix(x)
>             if(!is.numeric(z))
>                 stop("factor analysis applies only to numerical
variables")
>             if(!missing(subset)) z <- z[subset, , drop = FALSE]
>         }
>         covmat <- cov.wt(z)
>         cv <- covmat$cov
>         n.obs <- covmat$n.obs
>     }
>     else stop("'covmat' is of unknown type")
>     scores <- match.arg(scores)
>     if(scores != "none" && !have.x)
>         stop("requested scores without an 'x' matrix")
>     p <- ncol(cv)
>     if(p < 3) stop("factor analysis requires at least three variables")
>     dof <- 0.5 * ((p - factors)^2 - p - factors)
>     if(dof < 0)
>         stop(gettextf("%d factors is too many for %d variables", factors,
p),
>              domain = NA)
>     sds <- sqrt(diag(cv))
>     cv <- cv/(sds %o% sds)
> 
>     cn <- list(nstart = 1, trace = FALSE, lower = 0.005)
>     cn[names(control)] <- control
>     more <- list(...)[c("nstart", "trace", "lower", "opt", "rotate")]
>     if(length(more)) cn[names(more)] <- more
> 
>     if(is.null(start)) {
>         start <- (1 - 0.5*factors/p)/diag(solve(cv))
>         if((ns <- cn$nstart) > 1)
>             start <- cbind(start, matrix(runif(ns-1), p, ns-1,
byrow=TRUE))
>     }
>     start <- as.matrix(start)
>     if(nrow(start) != p)
>         stop(gettextf("'start' must have %d rows", p), domain = NA)
>     nc <- ncol(start)
>     if(nc < 1) stop("no starting values supplied")
> 
>     best <- Inf
>     for (i in 1:nc) {
>         nfit <- factanal.fit.mle(cv, factors, start[, i],
>                                  max(cn$lower, 0), cn$opt)
>         if(cn$trace)
>             cat("start", i, "value:", format(nfit$criteria[1]),
>                 "uniqs:", format(as.vector(round(nfit$uniquenesses, 4))),
> "\n")
>         if(nfit$converged && nfit$criteria[1] < best) {
>             fit <- nfit
>             best <- fit$criteria[1]
>         }
>     }
>     if(best == Inf) stop("unable to optimize from these starting
value(s)")
>     load <- fit$loadings
>     if(rotation != "none") {
>         rot <- do.call(rotation, c(list(load), cn$rotate))
>                                         # the following lines modified by
J.
> Fox, 26 June 2005
>         if (is.list(rot)){
>           load <- rot$loadings
>                                         #the following lines changed by
> Ulrich Keller, 9 Sept 2008
>           fit$rotmat <-
>             if(inherits(rot, "GPArotation")) {
>               t(solve(rot$Th))
>             } else {
>               rot$rotmat
>             }
>                                         #end changes Ulrich Keller, 9 Sept
> 2008
>         }
>         else load <- rot
>                                         # end modifications J. Fox, 26
June
> 2005
>     }
>     fit$loadings <- sortLoadings(load)
>     class(fit$loadings) <- "loadings"
>     fit$na.action <- na.act # not used currently
>     if(have.x && scores != "none") {
>         Lambda <- fit$loadings
>         zz <- scale(z, TRUE, TRUE)
>         switch(scores,
>                regression = {
>                    sc <- zz %*% solve(cv, Lambda)
>                    if(!is.null(Phi <- attr(Lambda, "covariance")))
>                        sc <- sc %*% Phi
>                },
>                Bartlett = {
>                    d <- 1/fit$uniquenesses
>                    tmp <- t(Lambda * d)
>                    sc <- t(solve(tmp %*% Lambda, tmp %*% t(zz)))
>                })
>         rownames(sc) <- rownames(z)
>         colnames(sc) <- colnames(Lambda)
>         if(!is.null(na.act)) sc <- napredict(na.act, sc)
>         fit$scores <- sc
>     }
>     if(!is.na(n.obs) && dof > 0) {
>         fit$STATISTIC <- (n.obs - 1 - (2 * p + 5)/6 -
>                      (2 * factors)/3) * fit$criteria["objective"]
>         fit$PVAL <- pchisq(fit$STATISTIC, dof, lower.tail = FALSE)
>     }
>     fit$n.obs <- n.obs
>     fit$call <- cl
>     fit
> }
> 
> factanal.fit.mle <-
>     function(cmat, factors, start=NULL, lower = 0.005, control = NULL,
...)
> {
>     FAout <- function(Psi, S, q)
>     {
>         sc <- diag(1/sqrt(Psi))
>         Sstar <- sc %*% S %*% sc
>         E <- eigen(Sstar, symmetric = TRUE)
>         L <- E$vectors[, 1:q, drop = FALSE]
>         load <- L %*% diag(sqrt(pmax(E$values[1:q] - 1, 0)), q)
>         diag(sqrt(Psi)) %*% load
>     }
>     FAfn <- function(Psi, S, q)
>     {
>         sc <- diag(1/sqrt(Psi))
>         Sstar <- sc %*% S %*% sc
>         E <- eigen(Sstar, symmetric = TRUE, only.values = TRUE)
>         e <- E$values[-(1:q)]
>         e <- sum(log(e) - e) - q + nrow(S)
> ##        print(round(c(Psi, -e), 5))  # for tracing
>         -e
>     }
>     FAgr <- function(Psi, S, q)
>     {
>         sc <- diag(1/sqrt(Psi))
>         Sstar <- sc %*% S %*% sc
>         E <- eigen(Sstar, symmetric = TRUE)
>         L <- E$vectors[, 1:q, drop = FALSE]
>         load <- L %*% diag(sqrt(pmax(E$values[1:q] - 1, 0)), q)
>         load <- diag(sqrt(Psi)) %*% load
>         g <- load %*% t(load) + diag(Psi) - S
>         diag(g)/Psi^2
>     }
>     p <- ncol(cmat)
>     if(is.null(start))
>         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),
>                  q = factors, S = cmat)
>     Lambda <- FAout(res$par, cmat, factors)
>     dimnames(Lambda) <- list(dimnames(cmat)[[1]],
>                              paste("Factor", 1:factors, sep = ""))
>     p <- ncol(cmat)
>     dof <- 0.5 * ((p - factors)^2 - p - factors)
>     un <- res$par
>     names(un) <- colnames(cmat)
>     class(Lambda) <- "loadings"
>     ans <- list(converged = res$convergence == 0,
>                 loadings = Lambda, uniquenesses = un,
>                 correlation = cmat,
>                 criteria = c(objective = res$value, counts = res$counts),
>                 factors = factors, dof = dof, method = "mle")
>     class(ans) <- "factanal"
>     ans
> }
> 
> print.loadings <- function(x, digits = 3, cutoff = 0.1, sort = FALSE, ...)
> {
>     Lambda <- unclass(x)
>     p <- nrow(Lambda)
>     factors <- ncol(Lambda)
>     if (sort) {
>         mx <- max.col(abs(Lambda))
>         ind <- cbind(1:p, mx)
>         mx[abs(Lambda[ind]) < 0.5] <- factors + 1
>         Lambda <- Lambda[order(mx, 1:p),]
>     }
>     cat("\nLoadings:\n")
>     fx <- format(round(Lambda, digits))
>     names(fx) <- NULL
>     nc <- nchar(fx[1], type="c")
>     fx[abs(Lambda) < cutoff] <- paste(rep(" ", nc), collapse = "")
>     print(fx, quote = FALSE, ...)
>     vx <- colSums(x^2)
>     varex <- rbind("SS loadings" = vx)
>     if(is.null(attr(x, "covariance"))) {
>         varex <- rbind(varex, "Proportion Var" = vx/p)
>         if(factors > 1)
>             varex <- rbind(varex, "Cumulative Var" = cumsum(vx/p))
>     }
>     cat("\n")
>     print(round(varex, digits))
>     invisible(x)
> }
> 
> 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, ...)
>                                         # 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) <- paste("Factor", 1:factors, sep="")
> 
>                                         # the following line changed by
> Ulrich Keller, 9 Sept 2008
>       if (TRUE != 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)) {
>         factors <- x$factors
>         cat("\nTest of the hypothesis that", factors, if(factors == 1)
>             "factor is" else "factors are", "sufficient.\n")
>         cat("The chi square statistic is", round(x$STATISTIC, 2), "on",
> x$dof,
>             if(x$dof == 1) "degree" else "degrees",
>             "of freedom.\nThe p-value is", signif(x$PVAL, 3), "\n")
>     } else {
>         cat(paste("\nThe degrees of freedom for the model is",
>                   x$dof, "and the fit was", round(x$criteria["objective"],
> 4),
>                   "\n"))
>     }
>     invisible(x)
> }
> 
> varimax <- function(x, normalize = TRUE, eps = 1e-5)
> {
>     nc <- ncol(x)
>     if(nc < 2) return(x)
>     if(normalize) {
>         sc <- sqrt(drop(apply(x, 1, function(x) sum(x^2))))
>         x <- x/sc
>     }
>     p <- nrow(x)
>     TT <- diag(nc)
>     d <- 0
>     for(i in 1:1000) {
>         z <- x %*% TT
>         B  <- t(x) %*% (z^3 - z %*% diag(drop(rep(1, p) %*% z^2))/p)
>         sB <- La.svd(B)
>         TT <- sB$u %*% sB$vt
>         dpast <- d
>         d <- sum(sB$d)
>         if(d < dpast * (1 + eps)) break
>     }
>     z <- x %*% TT
>     if(normalize) z <- z * sc
>     dimnames(z) <- dimnames(x)
>     class(z) <- "loadings"
>     list(loadings = z, rotmat = TT)
> }
> 
> promax <- function(x, m = 4)
> {
>     if(ncol(x) < 2) return(x)
>     dn <- dimnames(x)
>     xx <- varimax(x)
>     x <- xx$loadings
>     Q <- x * abs(x)^(m-1)
>     U <- lm.fit(x, Q)$coefficients
>     d <- diag(solve(t(U) %*% U))
>     U <- U %*% diag(sqrt(d))
>     dimnames(U) <- NULL
>     z <- x %*% U
>     U <- xx$rotmat %*% U
>     dimnames(z) <- dn
>     class(z) <- "loadings"
>     list(loadings = z, rotmat = U)
> }
> 
> --=-hiYzUeWcRJ/+kx41aPIZ
> Content-Disposition: attachment; filename="factanal.Rd"
> Content-Type: text/x-matlab; name="factanal.Rd"; charset="utf-8"
> Content-Transfer-Encoding: 8bit
> 
> % File src/library/stats/man/factanal.Rd
> % Part of the R package, http://www.R-project.org
> % Copyright 1995-2007 R Core Development Team
> % Distributed under GPL 2 or later
> 
> \name{factanal}
> \alias{factanal}
> %\alias{factanal.fit.mle}
> \encoding{latin1}
> \title{Factor Analysis}
> \description{
>   Perform maximum-likelihood factor analysis on a covariance matrix or
>   data matrix.
> }
> \usage{
> factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
>          subset, na.action, start = NULL,
>          scores = c("none", "regression", "Bartlett"),
>          rotation = "varimax", control = NULL, \dots)
> }
> \arguments{
>   \item{x}{A 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
>     \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
>     \code{\link{cov.wt}}.  Of course, correlation matrices are covariance
>     matrices.}
>   \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
>     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{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
>     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,
>     \describe{
>       \item{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
> 	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{\dots}{Components of \code{control} can also be supplied as
>     named arguments to \code{factanal}.}
> }
> \details{
>   The factor analysis model is
>   \deqn{x = \Lambda f + e}{ x = Lambda f + e}
>   for a \eqn{p}--element row-vector \eqn{x}, a \eqn{p \times k}{p x k}
>   matrix of \emph{loadings}, a \eqn{k}--element vector of \emph{scores}
>   and a \eqn{p}--element vector of errors.  None of the components
>   other than \eqn{x} is observed, but the major restriction is that the
>   scores be uncorrelated and of unit variance, and that the errors be
>   independent with variances \eqn{\Phi}{Phi}, the
>   \emph{uniquenesses}.  Thus factor analysis is in essence a model for
>   the covariance matrix of \eqn{x},
>   \deqn{\Sigma = \Lambda^\prime\Lambda + \Psi}{Sigma = Lambda'Lambda +
Psi}
>   There is still some indeterminacy in the model for it is unchanged
>   if \eqn{\Lambda}{Lambda} is replaced by \eqn{G\Lambda}{G Lambda} for
>   any orthogonal matrix \eqn{G}.  Such matrices \eqn{G} are known as
>   \emph{rotations} (although the term is applied also to non-orthogonal
>   invertible matrices).
> 
>   If \code{covmat} is supplied it is used.  Otherwise \code{x} is used if
>   it is a matrix, or a formula \code{x} is used with \code{data} to
>   construct a model matrix, and that is used to construct a covariance
>   matrix.  (It makes no sense for the formula to have a response,
>   and all the variables must be numeric.)  Once a covariance matrix is
found
> or
>   calculated from \code{x}, it is converted to a correlation matrix for
>   analysis.  The correlation matrix is returned as component
>   \code{correlation} of the result.
> 
>   The fit is done by optimizing the log likelihood assuming multivariate
>   normality over the uniquenesses.  (The maximizing loadings for given
>   uniquenesses can be found analytically: Lawley & Maxwell (1971,
>   p. 27).)  All the starting values supplied in \code{start} are tried
>   in turn and the best fit obtained is used.  If \code{start = NULL}
>   then the first fit is started at the value suggested by
>   \enc{J?reskog}{Joreskog} (1963) and given by Lawley & Maxwell
>   (1971, p. 31), and then \code{control$nstart - 1} other values are
>   tried, randomly selected as equal values of the uniquenesses.
> 
>   The uniquenesses are technically constrained to lie in \eqn{[0, 1]},
>   but near-zero values are problematical, and the optimization is
>   done with a lower bound of \code{control$lower}, default 0.005
>   (Lawley & Maxwell, 1971, p. 32).
> 
>   Scores can only be produced if a data matrix is supplied and used.
>   The first method is the regression method of Thomson (1951), the
>   second the weighted least squares method of Bartlett (1937, 8).
>   Both are estimates of the unobserved scores \eqn{f}.  Thomson's method
>   regresses (in the population) the unknown \eqn{f} on \eqn{x} to yield
>   \deqn{\hat f = \Lambda^\prime \Sigma^{-1} x}{hat f = Lambda' Sigma^-1 x}
>   and then substitutes the sample estimates of the quantities on the
>   right-hand side.  Bartlett's method minimizes the sum of squares of
>   standardized errors over the choice of \eqn{f}, given (the fitted)
>   \eqn{\Lambda}{Lambda}.
> 
>   If \code{x} is a formula then the standard NA-handling is applied to
>   the scores (if requested): see \code{\link{napredict}}.
> }
> \value{
>   An object of class \code{"factanal"} with components
>   \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.}
>   \item{uniquenesses}{The uniquenesses computed.}
>   \item{correlation}{The correlation matrix used.}
>   \item{criteria}{The results of the optimization: the value of the
>     negative log-likelihood and information on the iterations used.}
>   \item{factors}{The argument \code{factors}.}
>   \item{dof}{The number of degrees of freedom of the factor analysis
model.}
>   \item{method}{The method: always \code{"mle"}.}
>   %Modification by Ulrich Keller, Sept 9 2008
>   \item{rotmat}{The rotation matrix if relevant.}
>   %End modicication by Ulrich Keller, Sept 9 2008
>   \item{scores}{If requested, a matrix of scores.  \code{napredict} is
>     applied to handle the treatment of values omitted by the
> \code{na.action}.}
>   \item{n.obs}{The number of observations if available, or \code{NA}.}
>   \item{call}{The matched call.}
>   \item{na.action}{If relevant.}
>   \item{STATISTIC, PVAL}{The significance-test statistic and P value, if
>     if can be computed.}
> }
> 
> \note{
>   There are so many variations on factor analysis that it is hard to
>   compare output from different programs.  Further, the optimization in
>   maximum likelihood factor analysis is hard, and many other examples we
>   compared had less good fits than produced by this function.  In
>   particular, solutions which are Heywood cases (with one or more
>   uniquenesses essentially zero) are much often common than most texts
>   and some other programs would lead one to believe.
> }
> 
> \references{
>   Bartlett, M. S. (1937) The statistical conception of mental factors.
>   \emph{British Journal of Psychology}, \bold{28}, 97--104.
> 
>   Bartlett, M. S. (1938) Methods of estimating mental
>   factors. \emph{Nature}, \bold{141}, 609--610.
> 
>   \enc{J?reskog}{Joreskog}, K. G. (1963)
>   \emph{Statistical Estimation in Factor Analysis.}  Almqvist and
Wicksell.
> 
>   Lawley, D. N. and Maxwell, A. E. (1971) \emph{Factor Analysis as a
>     Statistical Method.} Second edition. Butterworths.
> 
>   Thomson, G. H. (1951) \emph{The Factorial Analysis of Human Ability.}
>   London University Press.
> }
> 
> \seealso{
>   \code{\link{print.loadings}},
>   \code{\link{varimax}}, \code{\link{princomp}},
>   \code{\link[datasets]{ability.cov}},
\code{\link[datasets]{Harman23.cor}},
>   \code{\link[datasets]{Harman74.cor}}
> }
> 
> \examples{
> # A little demonstration, v2 is just v1 with noise,
> # and same for v4 vs. v3 and v6 vs. v5
> # Last four cases are there to add noise
> # and introduce a positive manifold (g factor)
> v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6)
> v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)
> v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6)
> v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)
> v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5)
> v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)
> m1 <- cbind(v1,v2,v3,v4,v5,v6)
> cor(m1)
> factanal(m1, factors=3) # varimax is the default
> factanal(m1, factors=3, rotation="promax")
> # The following shows the g factor as PC1
> prcomp(m1)
> 
> ## formula interface
> factanal(~v1+v2+v3+v4+v5+v6, factors = 3,
>          scores = "Bartlett")$scores
> 
> ## a realistic example from Barthlomew (1987, pp. 61-65)
> utils::example(ability.cov)
> }
> \keyword{multivariate}
> 
> --=-hiYzUeWcRJ/+kx41aPIZ--
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel



[Rd] Addendum to wishlist bug report #10931 (factanal) (PR#12754)
John Fox