[R] bug in promax?

[Christof Schuster]

I was wondering whether the following inconsistency of the promax
rotation function with the results of a promax rotation using SAS
should be considered a bug in the promax function of R. Any comments
will be highly appreciated.

The following is a loading matrix obtained from a varimax rotation in SAS:
# Factor loadings after varimax rotation
x <- t(array(c(0.78107,         0.35573,
               0.85995,         0.07887,
               0.79376,         0.22499,
               0.26307,         0.85087,
               0.29449,         0.72966,
               0.06064,         0.68135), c(2,6)))

Clearly, the varimax function of R confirms this because

varimax(x)$rotmat
              [,1]         [,2]
[1,]  1.000000e+00 2.021128e-07
[2,] -2.021128e-07 1.000000e+00

yields the identity matrix.

SAS will return the following matrix after a promax rotation (for which
m=3 by default):

 0.75760         0.18383
 0.91956        -0.14061
 0.80720         0.03711
 0.05546         0.86360
 0.12295         0.72239
-0.11977         0.73116

However, using the promax function of R yields

promax(x, m=3)$loadings
            [,1]        [,2]
[1,]  0.74582399  0.20809487
[2,]  0.90713821 -0.11197901
[3,]  0.79547453  0.06260245
[4,]  0.05020993  0.86729962
[5,]  0.11747045  0.72788738
[6,] -0.12182636  0.72904497

Although the values are quite close to the SAS solution, there is
nevertheless a small mismatch. The following slightly edited function
of promax

mypromax <- function (x, m = 3) 
{
    if (ncol(x) < 2) 
        return(x)
    dn <- dimnames(x)
    xx <- varimax(x)
    x <- xx$loadings
    x1 <- diag(1/sqrt(diag(x %*% t(x)))) %*% x
    Q <- x1 * abs(x1)^(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
    list(loadings = z, rotmat = U)
}

yields 

round(mypromax(x)$loadings, 5)
         [,1]     [,2]
[1,]  0.75760  0.18383
[2,]  0.91956 -0.14061
[3,]  0.80721  0.03711
[4,]  0.05546  0.86360
[5,]  0.12295  0.72238
[6,] -0.11977  0.73116

Clearly, these values are virtually identical to the loadings after
promax rotation from SAS. The only differences of the above mypromax
function to the original function is that the loadings after varimax
rotation are first normalized (see the line in which x1 is calculated)
before they are used to calculate the target matrix Q in the next
line.

Hendrickson & White (1964), BJSP, 17, 65-70 who originally suggested
PROMAX wrote about the target matrix, p. 66: "Each element of this
matrix [the target matrix Q] is, ..., the kth power of the
corresponding element in the row-column normalized orthogonal matrix."

Although I am not sure what Hendrickson and White meant with "row-column"
normalization it appears to me that the original promax function is
missing a normalization of the loadings after varimax rotation before the
target matrix is calculated. The normalization I have used (and that
appears to be used by SAS) ensures that diag(x1 x1')=I, where x1 denotes
the normalized loadings.

Which set of loadings are the correct promax loadings?

Christof Schuster

Christof Schuster
University of Notre Dame
Department of Psychology                       
103 Haggar Hall
Notre Dame, IN 46556

Tel: (574) 631-5473    email: cschuste at nd.edu
Fax: (574) 631-8883    www.nd.edu/~cschuste