Alpha Factor Analysis

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: ALPHA                                               */
 /*   TITLE: Alpha Factor Analysis                               */
 /* PRODUCT: IML                                                 */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: MATRIX  STATAPP SUGI6                               */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: Rick Wicklin                UPDATE: Sep 2013        */
 /*     REF:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/


proc iml;
/*                Alpha Factor Analysis                      */
/*  Ref: Kaiser et al., 1965 Psychometrika, pp. 12-13        */
/*  Input:  r = correlation matrix                           */
/*  Output: m = eigenvalues                                  */
/*          h = communalities                                */
/*          f = factor pattern                               */
start alpha(m, h, f, r);
   p = ncol(r);
   q = 0;
   h = 0;                                      /* initialize */
   h2 = I(p) - diag(1/vecdiag(inv(r)));/* smc=sqrd mult corr */
   do while(max(abs(h-h2))>.001); /* iterate until converges */
      h = h2;
      hi = diag(sqrt(1/vecdiag(h)));
      g = hi*(r-I(p))*hi + I(p);
      call eigen(m,e,g);         /* get eigenvalues and vecs */
      if q=0 then do;
         q = sum(m>1);                  /* number of factors */
         iq = 1:q;
      end;                                   /* index vector */
      mm = diag(sqrt(m[iq,]));           /* collapse eigvals */
      e = e[,iq] ;                       /* collapse eigvecs */
      h2 = h*diag((e*mm) [,##]);        /* new communalities */
   end;
   hi = sqrt(h);
   h = vecdiag(h2);               /* communalities as vector */
   f = hi*e*mm;                         /* resulting pattern */
finish;

/* Correlation Matrix from Harmon, Modern Factor Analysis, */
/* Second edition, page 124, "Eight Physical Variables"    */
nm = {Var1 Var2 Var3 Var4 Var5 Var6 Var7 Var8};
r ={ 1.00 .846 .805 .859 .473 .398 .301 .382 ,
     .846 1.00 .881 .826 .376 .326 .277 .415 ,
     .805 .881 1.00 .801 .380 .319 .237 .345 ,
     .859 .826 .801 1.00 .436 .329 .327 .365 ,
     .473 .376 .380 .436 1.00 .762 .730 .629 ,
     .398 .326 .319 .329 .762 1.00 .583 .577 ,
     .301 .277 .237 .327 .730 .583 1.00 .539 ,
     .382 .415 .345 .365 .629 .577 .539 1.00};
run alpha(Eigenvalues, Communalities, Factors, r);
print Eigenvalues,
      Communalities[rowname=nm],
      Factors[label="Factor Pattern" rowname=nm];