;---------------------------------------------------------------------------- ; $Id: promax_hw64.pro,v 1.10 2002/11/06 02:59:39 johnny Exp $ ;+ ; NAME: ; PROMAX_HW64 ; ; PURPOSE: ; Given an input array of factor loadings (empirical orthogonal fun- ; ctions, i.e. EOFs), rotate the factors using a method almost iden- ; tical to the Promax method of Hendrickson and White (1964). ; ; CATEGORY: ; Statistics. ; ; CALLING SEQUENCE: ; Result = PROMAX_HW64(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. ; One can generate in_a by taking the transpose of the ; values returned in the Coefficients keyword in the built- ; in function PCOMP. Floating array dimensioned (r,p). ; Unchanged by function. ; ; INPUT KEYWORD PARAMETERS: ; M: Power m to raise "target matrix" to (while retaining sign ; of said target matrix) to create the "ideal pattern ; matrix" (RJ93) which decreases high loadings less than ; low loadings. See eqn. 28 in R86. If keyword is not ; defined, m defaults to 2 (following R86, p. 320). Un- ; changed by function. ; ; OUTPUT KEYWORD PARAMETERS: ; P_P: Column vectors of P_P are the first r Promax oblique ; primary factor pattern EOFs (factor loadings), floating ; array dimensioned (r,p). This is P_p in RJ93, p. 221. ; Created/overwritten by function. ; ; P_R: Column vectors of P_R are the first r Promax oblique ; reference pattern EOFs (factor loadings), floating array ; dimensioned (r,p). This is P_r in RJ93, p. 221. Created/ ; overwritten by function. ; ; PHI_P: Correlation between primary factors, floating array di- ; mentioned (r,r). Each [i,j] element is the correlation ; between factors i and j. This is phi_p in RJ93, p. 221. ; Created/overwritten by function. ; ; PHI_R: Correlation between reference axes, floating array dimen- ; tioned (r,r). Each [i,j] element is the correlation ; between reference axes i and j. This is phi_r in RJ93, ; p. 221. Created/overwritten by function. ; ; S_P: Column vectors of S_P are the first r Promax oblique ; primary factor structure EOFs (factor loadings), floating ; array dimensioned (r,p). This is S_p in RJ93, p. 221. ; Created/overwritten by function. ; ; S_R: Column vectors of S_R are the first r Promax oblique ; reference structure EOFs (factor loadings), floating ; array dimensioned (r,p). This is S_r in RJ93, p. 220 ; and is the same as the matrix passed out via Result. ; Created/overwritten by function. ; ; T_P: Primary structure transformation matrix, floating ; array dimensioned (r,r). This is T_p in RJ93, p. 221. ; Created/overwritten by function. ; ; T_R: Reference structure transformation matrix, floating ; array dimensioned (r,r). This is T_r in RJ93, p. 220. ; Created/overwritten by function. ; ; OUTPUTS: ; Result: Column vectors of Result are the first r Promax oblique ; reference structure EOFs (factor loadings), floating ; array dimensioned (r,p). This is S_r in RJ93, p. 220 ; and is the same as the matrix passed out by output ; keyword S_R. Created/overwritten by function. ; ; FILE DEVICE I/O: ; None. ; ; COMMON BLOCKS: ; None. ; ; EXAMPLE: ; This example uses input data from Table 7.I in RJ93 and gives ; results that are similar to Tables 7.IV and 7.V. The differences ; in output with those tables are probably because the Varimax routine ; used in this function probably has different iteration criteria than ; RJ93, and it appears in calculating Tables 7.IV and 7.V the B matrix ; was not forced row normalized in RJ93 (pp. 220, 343, and see note ; below in this header block). ; ; Given input factors, conduct a Promax 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] ] ; S_r = PROMAX_HW64(A, M=4, P_P=P_p, S_P=S_p) ; IDL returns for P_p: ; 0.0108750 0.912841 -0.160862 ; 1.02958 0.162965 0.254964 ; -0.0149609 0.168635 -0.910244 ; 0.902053 -0.182991 -0.301612 ; For S_p: ; -0.00528595 0.989576 -0.601337 ; 0.975345 -0.0105829 -0.00968674 ; 0.141439 0.606637 -0.988544 ; 0.965725 -0.0828474 -0.376990 ; ; REFERENCES: ; - Hendrickson, A. E., and P. O. White (1964), "Promax: A quick ; method for rotation to oblique simple structure," British J. ; Statist. Psychol., Vol. 17, pp. 65-70. Abbreviated HW64 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. ; - Richman, M. B. (1986), "Rotation of principal components," Journal ; of Climatology, Vol. 6, pp. 293-335. Abbreviated R86 in this ; function. ; ; MODIFICATION HISTORY: ; - 18 Oct 2002: Orig. ver. Johnny Lin, CIRES/University of Colorado. ; Email: air_jlin@yahoo.com. Passed passably adequate tests. ; - 4 Nov 2002: Fix variance conservation problem. Passed visual ; and reasonable values tests. ; ; NOTES: ; - Written for IDL 5.3. ; - The output structure and pattern matrices are all based upon a ; row normalized target matrix B and column normalized reference ; and pattern transformation matrices. Thus, even if the input ; EOFs are non-normalized, the output EOFs will be "normalized" ; (see references for this function). ; - The form of the HW64's Promax method used is that as described in ; RJ93 and R86. The "target matrix" is chosen to be the Varimax ro- ; tated EOFs (factor loadings) matrix. ; - In R86, p. 320, the columns and rows of B (called A* in R86) are ; normalized. In testing RPCA_LA00, which uses PROMAX_HW64, it was ; found that normalizing columns of B causes total variance to not be ; conserved during rotation (whether one uses the method of R86 or ; regular normalization). Thus, in this implementation of Promax, ; only rows are normalized in B; this is the major deviation from ; the canonical Promax method. ; - RJ93, pp. 215-220 has a detailed discussion of what "primary" and ; "reference" factors are as well as "structure" and "pattern" loadings. ; - This function was tested using RJ93, Table 7.III as values for the ; Varimax target matrix (instead of calculating the Varimax matrix using ; VARIMAX_K58), and without row-column normalizing B at the beginning ; of the Promax procedure. This yielded results nearly identical to ; Tables 7.IV and 7.V, which is the basis of the explanation of the ; example discrepancy with RJ93 discussed above. ; - All keyword parameters are optional unless otherwise stated. ; - No procedures called with _Extra keyword are invoked. ; - Non-built-in functions and procedures called: VARIMAX_K58. ;- ; Copyright (c) 2002 Johnny Lin. For licensing, distribution conditions, ; and contact information, see http://www.johnny-lin.com/lib.html. ;---------------------------------------------------------------------------- FUNCTION PROMAX_HW64, in_A $ , M = in_m $ , P_P = out_P_p $ , P_R = out_P_r $ , PHI_P = out_phi_p $ , PHI_R = out_phi_r $ , S_P = out_S_p $ , S_R = out_S_r $ , T_P = out_T_p $ , T_R = out_T_r $ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- COMPILE_OPT IDL2 ON_ERROR, 0 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) in_type = dims[3] ;- set input data type if (dims[0] ne 2) then MESSAGE, 'error--bad input dimensions' if (ABS(r) gt p) then MESSAGE, 'error--bad r' if (N_ELEMENTS(in_m) eq 0) then begin ;- set value of exponent m (see m = 2 ; descrip. of keyword M; default endif else begin ; for m is 2 m = in_m endelse ; ---------- Do Varimax Rotation To Get Preliminary Promax Solution --------- ; ; Section input: A (undefined on section exit) ; Section output: B, H_B ; ; Selected variable key: ; A First r unrotated EOFs (factor loadings), dimensioned (r,p). ; EOFs (factor loadings) of orig. 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_B The square root of communality as an array in which each ; element corresponds to the elements of B. Dimensioned ; (r,p). B = VARIMAX_K58(TEMPORARY(A), H=H) H_B = REBIN(REFORM(H, 1,p), r,p) ; -------------------------- Find Promax Solution --------------------------- ; ; Section input: B, H_B ; Section output: B (changed), P_p, P_r, S_p, S_r, T_p, T_r ; ; The obliquely rotated EOFs are column vectors of the output (P_p, P_r, ; etc.) ; ; Selected variable key: ; B Target matrix (RJ93, p. 220) which are the Varimax rotated ; EOFs, dimensioned (r,p). On section input, the first r ; rotated EOFs (factor loadings) are column vectors of B and ; called A* in R86, p. 320. Rows forced to be normalized in ; this section. ; B_star "Ideal pattern matrix" (RJ93, p. 220), dimensioned (r,p). ; Calculating using method described in R86, p. 320 (R86 ; calls B_star "B"). Note that RJ93 (p. 220) incorrectly ; implies B_star=B^m, without retaining the sign of B (R86 ; is correct, if one compares back to HW64). ; inv_stat Flag storing whether any of the INVERT functions in the ; section did not work correctly. If non-zero then fctn. ; will stop executing. ; max_col Maximum of abs. abs. value in each col. of B (row normal- ; ized). Dimensioned (r,p) where each element corresponds ; one-to-one with B. ; P_p See description of keyword P_P. ; P_r See description of keyword P_R. ; S_p See description of keyword S_P. ; S_r See description of keyword S_R. ; T_p See description of keyword T_P. ; T_p_trans Transpose of T_p. ; T_r See description of keyword T_R. ; T_r_trans Transpose of T_r. ; ; Algorithm: Most notation is similar to RJ93, p. 220. Algorithms are ; from both RJ93 and R86; see the sub-sections for details. See the ; "Notes" section in the top header for additional information re the ; Promax terminology. ;- set initial inv_stat to 0: inv_stat = 0 ;- normalization of rows of B by H_B (B is called A* in R86): B = TEMPORARY(B) / H_B ;- calculate "ideal pattern matrix" (B_star in RJ93, p. 220; algorithm ; follows R86 eqn. 28): neg_pts = WHERE(B lt 0.0, cnt_neg) ;+ indices of which elements are pos_pts = WHERE(B ge 0.0, cnt_pos) ; what sign in the B matrix B_star = B^m if (cnt_neg gt 0) then B_star[neg_pts] = (ABS(B_star))[neg_pts] * (-1.0) if (cnt_pos gt 0) then B_star[pos_pts] = (ABS(B_star))[pos_pts] ;- calculate transformation matrix T_r (RJ93, p. 220) and normalize col. ; (RJ93, p. 343); also calculate T_r_trans: T_r = INVERT(TRANSPOSE(B) ## B, status) ## TRANSPOSE(B) ## B_star inv_stat = inv_stat > status tmp = MAKE_ARRAY(Dimension=SIZE(T_r, /Dimensions), Type=in_type) idx = LINDGEN((SIZE(T_r, /Dimensions))[0]) tmp[idx,idx] = (TRANSPOSE(T_r) ## T_r)[idx,idx] T_r = TEMPORARY(T_r) ## SQRT( INVERT(TEMPORARY(tmp), status) ) inv_stat = inv_stat > status T_r_trans = TRANSPOSE(T_r) ;- calculate transformation matrix T_p (RJ93, p. 220) with the rows of ; transpose of T_p (T_p_trans) normalized (RJ93, p. 343): T_p_trans = INVERT(T_r, status) inv_stat = inv_stat > status tmp = MAKE_ARRAY(Dimension=SIZE(T_p_trans, /Dimensions), Type=in_type) idx = LINDGEN((SIZE(T_p_trans, /Dimensions))[0]) tmp[idx,idx] = (T_p_trans ## TRANSPOSE(T_p_trans))[idx,idx] T_p_trans = INVERT(SQRT(TEMPORARY(tmp)), status) ## TEMPORARY(T_p_trans) inv_stat = inv_stat > status T_p = TRANSPOSE(T_p_trans) ;- calculate structures (RJ93, p. 221): primary factor structure S_p ; and oblique reference structure S_r, both dimensioned (r,p): S_p = B ## T_p S_r = B ## T_r ;- calculate patterns (RJ93, p. 221): primary factor pattern P_p and ; oblique reference pattern P_r, both dimensioned (r,p): P_p = B ## INVERT(T_p_trans, status) inv_stat = inv_stat > status P_r = B ## INVERT(T_r_trans, status) inv_stat = inv_stat > status ;- check if any of the matrix inverses are bad: if (inv_stat ne 0) then MESSAGE, 'error--bad inversions' ; ---------- Calc. Correl. Betw. Primary Factors and Reference Axes --------- ; ; Section input: T_p, T_p_trans, T_r, T_r_trans ; Section output: phi_p, phi_r ; ; Selected variable key: ; phi_p Correlations between primary factors, dimensioned (r,r). ; phi_r Correlations between reference axes, dimensioned (r,r). ; ; Algorithm: See RJ93, p. 221 for details. phi_p = T_p_trans ## T_p phi_r = T_r_trans ## T_r ; ------------------------------ Prepare Output ----------------------------- out_P_p = TEMPORARY(P_p) ;- keyword P_P (optional) out_P_r = TEMPORARY(P_r) ;- keyword P_R (optional) out_phi_p = TEMPORARY(phi_p) ;- keyword PHI_P (optional) out_phi_r = TEMPORARY(phi_r) ;- keyword PHI_R (optional) out_S_p = TEMPORARY(S_p) ;- keyword S_P (optional) out_S_r = TEMPORARY(S_r) ;- keyword S_R (optional) out_T_p = TEMPORARY(T_p) ;- keyword T_P (optional) out_T_r = TEMPORARY(T_r) ;- keyword T_R (optional) Result = out_S_r ;- main function output RETURN, Result END ;=== end of function === ; ========== end file ==========