STATISTICA Basic Program: FactorAn.stb

Products Services Solutions Support Downloads Customers Partners Company
Home I Electronic Statistics Textbook I Contact Us I Request Quote I Search
STATISTICA BASIC Program FactorAn.stb
{ This program allows you to perform a principal components analysis, based on correlations, covariances, or the moment matrix (also, very large problems can be analyzed, e.g., >300 variables). The program will also perform a varimax rotation of the factor space, and compute the factor scores.
Program written, modified, or edited at StatSoft, Inc.}
randomaccess;
NoDataFileVariableNames;
	redim var(nvars);
	delta:=.00001; {minimum rotational angle}
	maxiter:=50;   {maximum number of iterations for rotation}
	markload:=.5;  {cutoff for highlighting loadings}
	if SelectVariables1 ('Select Variables for the Analysis',
	  3, nvars, var, n, 'Variables:')=0 then stop;

{check for missing data}
	for i:=1 to n do begin
	 nvalid:=ValCount (data(1,var(i)), 1, ncases);
	 if nvalid<>ncases then begin
	  DisplayMessageBox (MB_ICONSTOP, 'Missing Data',
	   'Some values in the selected variables have missing data; use the data
management tools to remove or replace the cases with missing data.');
	  stop;
	 end;
	end;

	iret:=DisplayButtonBox ('Print Results?',
	  'Display only|Print to printer|Print to text window');
	if iret=0 then stop;
	iprint:=iret-1;

	imethod:=DisplayButtonBox ('Matrix to Analyze',
	  'Covariances|Raw Crossproducts|Correlations');
	if imethod=0 then stop;
{allocate memory}

	redim covar(n,n);
	redim ddata(ncases,n);
	redim work1(ncases);

{copy data to ddata}
	for i:=1 to n do begin
	  MatrixGetColumn (data, var(i), work1);
	  MatrixSetColumn (ddata, i, work1);
	end;

{compute matrix to analyze}
	if imethod=3 then begin
	 line01$:='Correlations';
	 MatrixCorrelations (ddata, 1, covar);
	end;
	if imethod=2 then begin
	 line01$:='Raw Crossproducts';
	 MatrixCrossProductOfDev (ddata, 0, covar);
	 MatrixElemDivide (covar, ncases, covar);
	end;
	if imethod=1 then begin
	 line01$:='Covariances';
	 MatrixCrossProductOfDev (ddata, 1, covar);
	 MatrixElemDivide (covar, ncases-1, covar);
	end;
{release memory for ddata}
	redim ddata(1,1);
{display matrx to analyze}
	handle:=NewScrollsheet (n, n, covar, line01$,
	  ?RowNames$, ?ColumnNames$);
	ScrollsheetSetColumnWidth (handle, 8, 1);
	ScrollsheetSetRowNameWidth (handle, 8);
	for i:=1 to n do begin
	 ScrollsheetSetColumnName (handle, i, ?Name1$, varname(var(i)));
	 ScrollsheetSetRowName (handle, i, varname(var(i)));
	end;
      if iprint=1 then PrintScrollsheet (handle)
	else if iprint=2 then PrintScrollsheetToOutputWindow (handle);

{perform factor analysis}
	redim evals(n);
	redim evectors(n,n);
	MatrixEigenvectors (covar, evals, evectors, neigen);

{display eigenvalues}
	redim table(neigen,4);
	ValSum (evals, 1, n, sumeigen);
	cumul:=0;
	cumulp:=0;
	for i:=1 to neigen do begin
	 table(i,1):=evals(i);
	 table(i,2):=100*evals(i)/sumeigen;
	 cumul:=cumul+evals(i);
	 cumulp:=cumulp+table(i,2);
	 table(i,3):=cumul;
	 table(i,4):=cumulp;
	end;
	handle:=NewScrollsheet (neigen, 4, table, 'Eigenvalues',
	 ?RowNames$, ?ColumnNames$);
	ScrollsheetSetColumnWidth (handle, 8, 1);
	ScrollsheetSetColumnName (handle, 1, ?Name1$, 'Eigenv.');
	ScrollsheetSetColumnName (handle, 2, 'Percent', 'Eigenv.');
	ScrollsheetSetColumnName (handle, 3, 'Cumulatv', 'Eigenv.');
	ScrollsheetSetColumnName (handle, 4, '% Cumul.', 'Eigenv.');
      if iprint=1 then PrintScrollsheet (handle)
	else if iprint=2 then PrintScrollsheetToOutputWindow (handle);

{plot of eigenvalues}
	redim work1(n);
	for i:=1 to n do work1(i):=i;
	handle:=NewGraph (lineplot, 'Plot of Eigenvalues',
	  'Eigenvalue', '', neigen, work1, evals);
	GraphSetPlotPointStyle (handle, 1, ON, ?Pattern, ?Size, ?Color);
	DisplayGraph (handle);
      if iprint=1 then PrintGraph(handle)
	else if iprint=2 then PrintGraphToOutputWindow(handle);


{retain how many factors?}
	nfactors:=2;
	DisplayNumericInputBox ('Retain how many factors?',
	 'Number of factors:', nfactors);

{compute unrotated factor loadings}
	redim loadings(n,n);
	redim table(n,n);
	redim work1(n);
	MatrixSetDiagonal (loadings, evals);
	for i:=1 to n do loadings(i,i):=sqrt(abs(loadings(i,i)));
	MatrixMultiply (evectors, loadings, loadings);

{check for signs; multiple by -1 so that most loadings are positive}
	MatrixIsGreater (loadings, table, table);
	for i:=1 to n do begin
	 MatrixGetColumn (table, i, work1);
	 ValSum (work1, 1, n, j);
	 if j<n/2 then begin
	  MatrixGetColumn (loadings, i, work1);
	  MatrixElemMultiply (work1, -1, work1);
	  MatrixSetColumn (loadings, i, work1);
	 end;
	end;

{display unrotated factor loadings}
	handle:=NewScrollsheet (n, nfactors, Loadings, 'Unrotated Factor Loadings',
	 ?RowNames$, ?ColumnNames$);
	ScrollsheetSetColumnWidth (handle, 8, 1);
	ScrollsheetSetRowNameWidth (handle, 8);
	for i:=1 to n do begin
	 if i<=nfactors then
	   ScrollsheetSetColumnName (handle, i, 'Factor','  '+str(i,4,0));
	 ScrollsheetSetRowName (handle, i, varname(var(i)));
	 for j:=1 to nfactors do begin
	  if abs(loadings(i,j))>markload then
	    ScrollsheetSetHilite (handle, i, j, 1);
	 end;
	end;
      if iprint=1 then PrintScrollsheet (handle)
	else if iprint=2 then PrintScrollsheetToOutputWindow (handle);

	if nfactors>1 then begin
{compute varimax (normalized) rotated factor solution}
{allocate memory for necessary arrays}
	for iter:=1 to maxiter do begin
	 redim u(n),vv(n),h(n),x(n),y(n);
	 redim l(n,2);
	 redim t(2,2);
	 MatrixSetToZero (h);
{square root of communalities}
	 for i:=1 to n do begin
	  for j:=1 to nfactors do h(i):=h(i)+loadings(i,j)^2;
	  h(i):=sqrt(h(i));
	 end;
{now perform varimax rotation of factor loadings}
{cycle through all pairs of factors}
	 rotated:=0;
	 for ifactor:=1 to nfactors do begin
	  for jfactor:=ifactor+1 to nfactors do begin
	   MatrixSetToZero (u);
	   MatrixSetToZero (vv);
	   MatrixSetToZero (x);
	   MatrixSetToZero (y);
	   MatrixSetToZero (l);
	   MatrixSetToZero (t);
	   c:=0;d:=0;
{normalize loadings for pair}
	   for i:=1 to n do begin
	    x(i):=loadings(i,ifactor)/h(i);
	    y(i):=loadings(i,jfactor)/h(i);
	   end;
	   for i:=1 to n do begin
	    u(i):=x(i)^2-y(i)^2;
	    vv(i):=2*x(i)*y(i);
	    c:=c+u(i)^2-vv(i)^2;
	    d:=d+u(i)*vv(i);
	   end;
	   d:=d*2;
	   ValSum (u,1,n,a);
	   ValSum (vv,1,n,b);
	   e:=d-2*a*b/n;
	   f:=c-(a^2-b^2)/n;
	   t4:=e/f;
	   xx:=ArcTan(t4);
	   if f<0 then begin
	    if e>=0 then xx:=xx+2*pi;
	    xx:=xx-pi;
	   end;
	   xx:=xx/4;
	   if abs(xx)>delta then begin
{rotate pair; compute transformation matrix}
	    rotated:=1;
	    t(1,1):=cos(xx);
	    t(1,2):=-sin(xx);
	    t(2,1):=sin(xx);
	    t(2,2):=cos(xx);
	    MatrixSetColumn (l, 1, x);
	    MatrixSetColumn (l, 2, y);
	    MatrixMultiply (l, t, l);
	    for i:=1 to n do begin
{undo standardization}
	     for j:=1 to 2 do l(i,j):=l(i,j)*h(i);
	    end;
	    MatrixCopy (l, 1, 1, n, 1, loadings, 1, ifactor);
	    MatrixCopy (l, 1, 2, n, 1, loadings, 1, jfactor);
	   end; {end rotation of pair}
	  end; {end jfactor}
	 end;  {end ifactor}
{stop iterations if no rotation occurred in this iteration}
	 if rotated=0 then iter:=maxiter;
	end; {next iteration}
{check for signs; multiple by -1 so that most loadings are positive}
	MatrixSetToZero (table);
	MatrixIsGreater (loadings, table, table);
	for i:=1 to nfactors do begin
	 MatrixGetColumn (table, i, work1);
	 ValSum (work1, 1, n, j);
	 if j<n/2 then begin
	  MatrixGetColumn (loadings, i, work1);
	  MatrixElemMultiply (work1, -1, work1);
	  MatrixSetColumn (loadings, i, work1);
	 end;
	end;

{display varimax rotated factor loadings}
	handle:=NewScrollsheet (n, nfactors, Loadings, 'Varimax Rotated Factor Loadings',
	 ?RowNames$, ?ColumnNames$);
	ScrollsheetSetColumnWidth (handle, 8, 1);
	ScrollsheetSetRowNameWidth (handle, 8);
	for i:=1 to n do begin
	 if i<=nfactors then
	   ScrollsheetSetColumnName (handle, i, 'Factor','  '+str(i,4,0));
	 ScrollsheetSetRowName (handle, i, varname(var(i)));
	 for j:=1 to nfactors do begin
	  if abs(loadings(i,j))>markload then
	    ScrollsheetSetHilite (handle, i, j, 1);
	 end;
	end;
      if iprint=1 then PrintScrollsheet (handle)
	else if iprint=2 then PrintScrollsheetToOutputWindow (handle);

{plot of factor loadings (optional)}
	 iret:=DisplayMessageBox (MB_YESNOCANCEL, 'Plot Factor Loadings?',
	 'Do you want to plot the factor loadings for all pairs of factors in
scatterplots?');
	 if iret=IDCANCEL then stop;
	 if iret=IDYES then begin
	  for ifactor:=1 to nfactors do begin
	   for jfactor:=ifactor+1 to nfactors do begin
	    MatrixExtract (loadings, 1, ifactor, n, 1, x);
	    MatrixExtract (loadings, 1, jfactor, n, 1, y);
	    handle:=NewGraph (SCATTERPLOT,
	     'Rotated Factor Loadings|Factor '+Str (ifactor, 2, 0)+
		' vs. Factor '+Str (jfactor, 2, 0),
            'Factor '+Str(jfactor,2,0),
		'Factor '+Str(ifactor,2,0),
		n, x, y);
	    GraphSetPlot2DLayout (handle, 1, SCATTERPLOT,
		DATALABELS_TEXT,
		?BarStyle, ?BarWidth, ?DevLevel, ?IsRightAxis);
	    if imethod=3 then begin
	     GraphSetScaling (handle, AX_X, SCALING_MANUAL_0,
			-1, +1, .2);
	     GraphSetScaling (handle, AX_Y, SCALING_MANUAL_0,
			-1, +1, .2);
	     x(1):=-1;x(2):=+1;x(3):=missing;x(4):=0;x(5):=0;
	     y(1):=0;y(2):=0;y(3):=missing;y(4):=-1;y(5):=+1;
	     GraphAddPlot (handle, LINEPLOT,
		'', 5, x, y);
	     GraphSetPlotLineStyle (handle, 2, ON, L_SOLID, ?Size, BLACK);
	    end;
	    DisplayGraph (handle);
          if iprint=1 then PrintGraph(handle)
	    else if iprint=2 then PrintGraphToOutputWindow(handle);
  	   end;
	  end;
	 end;
	end;   {end of rotation (if n factors>1)}

	if imethod<>3 then stop;
{compute factor score coefficients}
	redim table(n,nfactors);
	MatrixExtract (loadings, 1, 1, n, nfactors, table);
	redim loadings(n,nfactors);
	MatrixDuplicate (table, loadings);
	redim table(nfactors,n);
	MatrixTranspose (loadings, table);
	redim work2(nfactors,nfactors);
	MatrixMultiply (table, loadings, work2);
	redim table(n,nfactors);
	MatrixInverse (work2, work2);
	MatrixMultiply (loadings, work2, table);


{display factor score coefficients}
	line01$:='Factor Score Coefficients';
	if nfactors>1 then line01$:=line01$+' (Rotated Solution)';
	handle:=NewScrollsheet (n, nfactors, table,
	 line01$,
	 ?RowNames$, ?ColumnNames$);
	ScrollsheetSetColumnWidth (handle, 8, 1);
	ScrollsheetSetRowNameWidth (handle, 8);
	for i:=1 to n do begin
	 if i<=nfactors then
	   ScrollsheetSetColumnName (handle, i, 'Factor','  '+str(i,4,0));
	 ScrollsheetSetRowName (handle, i, varname(var(i)));
	end;
      if iprint=1 then PrintScrollsheet (handle)
	else if iprint=2 then PrintScrollsheetToOutputWindow (handle);


{compute optional factor scores}
	iret:=DisplayMessageBox (MB_OKCANCEL, 'Compute Factor Scores?',
	'Go ahead and compute factor scores for the current factor solution?');
	if iret=IDCANCEL then stop;
	redim ddata(ncases,n);
	redim work1(ncases);
	for i:=1 to ncases do begin
	 for j:=1 to n do ddata(i,j):=data(i,var(j));

	end;

{standardize/subtract mean if necessary}
	if imethod<>2 then begin
	 for i:=1 to n do begin
	  MatrixExtract (ddata, 1, i, ncases, 1, work1);
	  ValMean (work1, 1, ncases, mean);
	  ValStDeviation (work1, 1, ncases, stdv);
	  MatrixElemSubtract (work1, mean, work1);
	  if imethod=3 then MatrixElemDivide (work1, stdv, work1);
	  MatrixSetColumn (ddata, i, work1);
	 end;
	end;
	redim table2(ncases,nfactors);
	MatrixMultiply (ddata, table, table2);
	line01$:='Factor Scores';
	if nfactors>1 then line01$:=line01$+' (Rotated Solution)';
	handle:=NewScrollsheet (ncases, nfactors, table2,
	 line01$,
	 ?RowNames$, ?ColumnNames$);
	ScrollsheetSetColumnWidth (handle, 8, 1);
	ScrollsheetSetRowNameWidth (handle, 8);
	for i:=1 to ncases do begin
	 if i<=nfactors then
	   ScrollsheetSetColumnName (handle, i, 'Factor','  '+str(i,4,0));
	 ScrollsheetSetRowName (handle, i, CaseName (i));
	end;
	if iprint=1 then PrintScrollsheet (handle)
	else if iprint=2 then PrintScrollsheetToOutputWindow (handle);
Back to List of Programs
[StatSoft]
2300 East 14th Street, Tulsa, OK 74104
Phone: (918) 749-1119; Fax: (918) 749-2217
e-mail: info@statsoft.com