;---------------------------------------------------------------------------- ; $Id: varimax_k58.pro,v 1.3 2002/10/11 21:27:10 johnny Exp johnny $ ;+ ; NAME: ; VARIMAX_K58 ; ; PURPOSE: ; Given an input array of factor loadings (empirical orthogonal fun- ; ctions, i.e. EOFs), rotate the factors using the Varimax criteria ; of Kaiser (1958, 1959). ; ; CATEGORY: ; Statistics. ; ; CALLING SEQUENCE: ; Result = VARIMAX_K58(in_A) ; ; INPUTS: ; in_A: Column vectors of in_A are the first r unrotated EOFs ; (factor loadings) of the input data (of n observations ; at p locations) from which these EOFs are derived. ; Floating array dimensioned (r,p). Unchanged by function. ; ; INPUT KEYWORD PARAMETERS: ; None. ; ; OUTPUT KEYWORD PARAMETERS: ; H: Square root of communality, floating array dimensioned ; (p), of each location (variable). Created/overwritten ; by function. ; ; FRAC_VAR: Fraction of variance accounted for by each r rotated ; EOF (factor), floating array dimensioned (r), assuming ; that the locations (variables) were standardized (see ; p. 203, 213, and other sections in RJ93 for details). ; Created/overwritten by function. ; ; OUTPUTS: ; Result: Column vectors of Result are the first r rotated EOFs ; (factor loadings), floating array dimensioned (r,p). ; Created/overwritten by function. ; ; FILE DEVICE I/O: ; None. ; ; COMMON BLOCKS: ; None. ; ; EXAMPLE: ; This example is taken from Table 7.I in RJ93. The small differences ; in output with Table 7.III are likely due to the use of a different ; criteria in this function for determining when the Varimax solution ; has converged. ; ; Given input factors, conduct a Varimax rotation: ; A = [ [ 0.4320, 0.8129, 0.3872] $ ; , [ 0.7950, -0.5416, 0.2565] $ ; , [ 0.5944, 0.7234, -0.3441] $ ; , [ 0.8945, -0.3921, -0.1863] ] ; B = VARIMAX_K58(A, H=H, Frac_Var=frac_var) ; IDL returns for B: ; -0.0204423 0.938674 -0.340334 ; 0.983662 0.0730206 0.134997 ; 0.0826106 0.435975 -0.893379 ; 0.939901 -0.0965213 -0.309596 ; For H^2: ; 0.997354 0.991148 0.995024 0.988580 ; For fraction of variance: ; 0.464562 0.271458 0.257007 ; ; REFERENCES: ; - Kaiser, H. F. (1958), "The varimax criterion for analytic rotation ; in factor analysis," Psychometrika, Vol. 23, pp. 187-200. Abbrev- ; iated K58 in this function. ; - Kaiser, H. F. (1959), "Computer program for varimax rotation in ; factor analysis," Educ. Psych. Meas., Vol. 19, No. 3, pp. 413-420. ; Abbreviated K59 in this function. ; - Reyment, Richard A., and K. G. Joreskog (1993), Applied Factor ; Analysis in the Natural Sciences, 2nd Edition (Cambridge: Cambridge ; University Press), 371 pp. Abbreviated RJ93 in this function. See ; also MATLAB code by L. F. Marcus in the Appendix of this book. ; ; MODIFICATION HISTORY: ; - 11 Oct 2002: Orig. ver. Johnny Lin, CIRES/University of Colorado. ; Email: air_jlin@yahoo.com. Passed passably adequate tests. ; ; NOTES: ; - Written for IDL 5.3. ; - All keyword parameters are optional unless otherwise stated. ; - No procedures called with _Extra keyword are invoked. ; - Only built-in functions and procedures are called. ;- ; Copyright (c) 2002 Johnny Lin. For licensing, distribution conditions, ; and contact information, see http://www.johnny-lin.com/lib.html. ;---------------------------------------------------------------------------- FUNCTION VARIMAX_K58, in_A $ , H = out_H $ , FRAC_VAR = out_frac_var $ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- COMPILE_OPT IDL2 ON_ERROR, 0 epsilon = (MACHAR()).eps ;- set machine precision tolerance A = in_A ;- protect input data dims = SIZE(A) ;- dimensions of input data r = dims[1] ;- set number of retained EOFs (factor loadings) ; (in RJ93 eqn. 7.2 etc., this is called k): p = dims[2] ;- set number of locations (in RJ93 these are called ; variables) if (dims[0] ne 2) then MESSAGE, 'error--bad input dimensions' if (ABS(r) gt p) then MESSAGE, 'error--bad r' ; --------------------------- Do Varimax Rotation --------------------------- ; ; Section input: A (undefined upon section output) ; Section output: B, H, H_B ; ; Selected variable key: ; A Unrotated EOFs, dimensioned (r,p). EOFs (factor loadings) ; of input data are column vectors of A. ; B Rotated EOFs, dimensioned (r,p). The first r rotated EOFs ; (factor loadings) are column vectors of B. ; H Square root of communality, dimensioned (p) ; H_B Same as H, but as an array in which each element corres- ; ponds to the elements of B and is dimensioned (r,p). ; ; Note: Some of the K59 local variables have the same names as variables ; in this function (e.g. A). To prevent a mix-up, all the K59 inter- ; mediate variables that are one letter long (and not capital X or Y) are ; named in this function as "duplicated" characters, e.g. xx is the same ; as x in K59, AA is the same as A in K59, etc. ; ; Algorithm: Rotation is done 2 EOFs at a time until there is a conver- ; gence of solution (i.e. a set of rotated EOFs B in which the Varimax ; criteria is maximized). Following the idea in K59, convergence of ; solution is assumed when we have tried to rotate through all r*(r-1)/2 ; pairs of EOFs in the matrix but all rotation attempts have failed as ; the angle for all pairs is essentially zero. Then the flag not_con- ; verged is set to 0 (false) and the rotation cycle while loop ends. ;--> some initial settings: B = TEMPORARY(A) ;- set initial B to A not_converged = 1 ;- init. flag to Varimax solution has not yet converged ;--> calculate square root of communality of each location (variable): ; this is the square root of the sum of the squares of each row of B ; (i.e. sum of squares of loadings for each p location) since the A ; factors (which B currently equals) are uncorrelated (see RJ93, eqns. ; 3.14, 3.1): H = SQRT( TOTAL(B^2, 1) ) H_B = REBIN(REFORM(H, 1,p), r,p) ;--> calc. B normalized by square root of communality (i.e. b_{ij}/h_{i}, ; e.g. in eqn. 7.3 in RJ93) and set as array bh, dimensioned (r,p): bh = B / H_B ;--> iterate until solution has converged (loop begin): while not_converged do begin ;- set/reset counter of number of r*(r-1)/2 pairs of EOFs (num_pairs): num_pairs = r * (r-1) / 2 ;- rotate all factors (EOFs) in B two factors at a time: each factor index ; i [0:r-1-1] is paired with each of the remaining factors (index j): for i=0,r-1-1 do begin ;+ go through factor index 0:r-1-1 (begin) for j=i+1,r-1 do begin ;+ pair i to "rest" of factors (begin) ;> select pair of factors at factors i,j (notation closely follows K59): xx = bh[i,*] yy = bh[j,*] ;> calculate angle of rotation (phi) that maximizes the Varimax function ; for factors i,j (notation closely follows K59; n in K59 is p here; ; see also p. 340 in RJ93); the algorithm used here to calculate phi ; already accounts for maximizing the Varimax function: uu = xx^2 - yy^2 vv = 2.0 * xx * yy AA = TOTAL(uu) BB = TOTAL(vv) CC = TOTAL(uu^2 - vv^2) DD = TOTAL(2.0 * uu * vv) num = DD - (2.0 * AA * BB / p) den = CC - ( (AA^2 - BB^2) / p ) phi = ATAN(num, den) / 4.0 ;> if phi is non-zero, rotate to get new factors (xx_rot and yy_rot, ; which are X and Y in K59), insert into bh; if phi is zero, decrement ; num_pairs by 1 and stop rotation process after one pass through all ; EOF pairs in bh have been shown to have phi=0 rotation angles: if (ABS(phi) ge 2.0*epsilon) then begin xx_rot = ( COS(phi) * xx) + (SIN(phi) * yy) yy_rot = (-SIN(phi) * xx) + (COS(phi) * yy) bh[i,*] = TEMPORARY(xx_rot) bh[j,*] = TEMPORARY(yy_rot) endif else begin num_pairs = num_pairs - 1 if (num_pairs eq 0) then not_converged = 0 endelse endfor ;+ pair i to "rest" of factors (end) endfor ;+ go through factor index 0:r-1-1 (end) ;--> iterate until solution has converged (loop end): endwhile ;--> denormalize (and deallocate) bh to get final B by multiplying bh by ; the square root of communality: B = TEMPORARY(bh) * H_B ; ---------------------- Calculate Fraction of Variance --------------------- ; ; Section input: B ; Section output: fvar ; ; Algorithm: Calculate fraction of variance (fvar) for each EOF/PC ; (factor) for the B matrix, as described on the bottom of p. 213, ; Table 7.I-7.III, and top of p. 215 of RJ93. fvar = TOTAL(B^2, 2) / FLOAT(p) ; ------------------------------ Prepare Output ----------------------------- out_H = TEMPORARY(H) ;- output sqrt of communality (optional) out_frac_var = TEMPORARY(fvar) ;- output frac. of var. of EOF/PCs (optional) Result = TEMPORARY(B) ;- output rotated EOFs RETURN, Result END ;=== end of function === ; ========== end file ==========