PRINFACTORS - updated version

[TomDoan]

This has been updated to include several new options. It also sign-flips loading vectors which are predominantly negative.

Code: Select all

*
*  PRINFACTORS does a principal components based factor analysis
*  of an input covariance or correlation matrix.
*
*  Syntax:
*
*  @PRINFACTORS( options ) sigma
*
*  Options:
*    NCOMPS=number of components [number of variables]
*    PRINT/[NOPRINT]  Controls whether the output will be displayed.
*    SIGMAHAT=output sigma matrix [not used]
*    LOADINGS=RECTANGULAR of loadings of components onto variables [not used]
*    VECTORS=RECTANGULAR of eigenvectors [not used]
*    VALUES=VECTOR of eigenvalues [not used]
*    COMMUNALITIES=VECTOR of communalities for variables [not used]
*
*  Parameters:
*    sigma -- The input covariance matrix.
*
*  Revision Schedule
*
*  05/05 Written by Tom Doan, Estima
*  08/06 Corrections to the comments describing the procedure syntax.
*  03/07 Define %VARIANCE. Add COMMUNALITIES, VECTORS and VALUES option.
*        Sign flip loading vectors which are primarily negative.
*
procedure prinfactors sigma
type symm sigma
*
option integer ncomps   0
option switch  print    0
option symm             *sigmahat
option rect             *loadings
option rect             *vectors
option vect             *values
option vect             *communalities
*
local rect eigenvec
local vect eigenval princoef scalings
local int  n nvars icomp
local real totalvar running
*
compute nvars=%rows(sigma)
*
* Do an eigen decomposition of the input matrix.
*
eigen sigma eigenval eigenvec
if %defined(values)
   compute values=eigenval
if %defined(vectors)
   compute vectors=eigenvec
*
* Figure out how many components need to be computed
*
if .not.ncomps
   compute n = nvars
else
   compute n = ncomps

if print {
   report(action=define)
   report(atrow=1,atcol=2,span) "Principal Components Analysis"
   report(atcol=1,atrow=2,fillby=cols) "Eigenvalue" "Proportion" "Cumulative" "Eigenvector"
}
*
* Get the total variation for normalizing the contributions
*
compute totalvar=%sum(eigenval)
*
* Initialize the running sum of the covariances
*
if %defined(sigmahat)
   compute sigmahat=%zeros(nvars,nvars)
*
* And the running sum of the variance explained
*
if %defined(communalities)
   compute communalities=%zeros(nvars,1)

compute running=0.0
do icomp=1,n
   compute running=running+eigenval(icomp)
   compute princoef=%xcol(eigenvec,icomp)
   if %defined(communalities)
      compute communalities=communalities+eigenval(icomp)*princoef.^2
   if print
      report(col=new,atrow=2,fillby=cols) eigenval(icomp) eigenval(icomp)/totalvar running/totalvar princoef
   if %defined(sigmahat)
      compute sigmahat=sigmahat+%outerxx(princoef)*eigenval(icomp)
end do icomp
if print
  report(action=show)
*
* Replace the diagonal elements of sigmahat with the input elements
*
if %defined(sigmahat) {
   do icomp=1,nvars
      compute sigmahat(icomp,icomp)=sigma(icomp,icomp)
   end do i
}
compute %variance   =totalvar
*
*  Compute the loadings if requested. Sign flip column which is predominantly negative
*
if %defined(loadings) {
   dim loadings(nvars,n)
   dim scalings(n)
   do j=1,n
      compute scalings(j)=sqrt(eigenval(j))*%sign(%sum(%xcol(eigenvec,j)))
   end do j
   ewise loadings(i,j)=eigenvec(i,j)*scalings(j)
}
*
end

[dniggeler@gmx.ch]

Hi

Could you please explain to a Rats beginner why in contrast to other stats program I get eigenvectors with a flipped sign.

Thank you for your help,

Dieter.

[TomDoan]

The signs of eigenvectors are arbitrary. RATS uses a translation of EISPACK FORTRAN code which is probably almost identical to the newer LAPACK code. There's no sign choice in that at all; the signs of the eigenvectors are what they are based upon the input matrix and the order of calculations.

The sign choice is mathematically irrelevant - it only affects how easy it is to read the results. What sign choice is most useful will actually depend upon the situation. The recoding of PRINFACTORS makes the sum of the elements in an eigenvector positive. While this is arguably more useful than having the sign determined effectively randomly, I could think of situations where making the largest value positive would make it easier to read. (CATS, for instance, by default will normalize an eigenvector to make its largest element 1.0).

[IRJ]

How can one extract the principal components from the procedure @prinfactors? More specifically, how can one obtain the output of the procedure @princomp using @prinfactors? And how can one normalize the eigenvectors in the two procedures so that the largest eigenvector takes a value of one?

[TomDoan]

How can one extract the principal components from the procedure @prinfactors? More specifically, how can one obtain the output of the procedure @princomp using @prinfactors? And how can one normalize the eigenvectors in the two procedures so that the largest eigenvector takes a value of one?

If you want the principal components, just use @PRINCOMP.

As an example of rescaling:

compute eigen=%ranmat(5,5)
dec vect maxv(%rows(eigen))
ewise maxv(i)=%maxvalue(%abs(%xcol(eigen,i)))
compute eigen=%ddivide(eigen,maxv)