1 1 MODEL SELECTION ALGORITHM FOR THE ORTHOGONAL FACTOR MODEL USING FACAIC HAMPARSUM BOZDOGAN AND DONALD E. RAMIREZ DEPARTMENT OF MATHEMATICS MATH-ASTRO BUILDING UNIVERSITY OF VIRGINIA CHARLOTTESVILLE, VIRGINIA 22903 TECHNICAL REPORT 16 (DECEMBER 3, 1986) 1 INTRODUCTION THIS TECHNICAL REPORT PRESENTS THE PROGRAM FACAIC WHICH IS DEVELOPED IN OUR PAPER: FACAIC: MODEL SELECTION ALGORITHM FOR THE ORTHOGONAL FACTOR MODEL USING AIC AND CAIC, PSYCHOMETRIKA (TO APPEAR) FACAIC IS A FORTRAN 77 PROGRAM WHICH FINDS THE BEST NUMBER OF FACTORS IN AN ORTHOGONAL FACTOR MODEL. THE CRITERIA ALLOWED ARE AKAIKE'S INFORMATION CRITERION AIC AND ITS CONSISTENT VARIATION CAIC. THE USER IS ASSUMED TO HAVE COMPUTED THE CORRELATION MATRIX R FROM THE DATA MATRIX X. THE USER INPUTS TO THE SUBROUTINE FACAIC ARE N THE NUMBER OF CASES, P THE NUMBER OF VARIABLES, DOAIC .TRUE. FOR AIC, CORR THE CORRELATION MATRIX IN SYMMETRIC STORAGE, AND SETNF .FALSE. FOR AN AUTOMATIC SEARCH. INDEX APPENDIX A CONTAINS THE RAW DATA SET WITH N=305 AND P=7 APPENDIX B CONTAINS THE TERMINAL OUTPUT FROM FACAIC APPENDIX C CONTAINS THE FULL OUTPUT FROM THE PROGRAM FACAIC APPENDIX D CONTAINS THE FORTRAN PROGRAM FACAIC APPENDIX E CONTAINS THE CALL STATEMENTS FROM FACAIC ************************************************************************ APPENDIX A: THE ORIGINAL DATA 1.0 .805 1.0 .859 .801 1.0 .473 .380 .436 1.0 .398 .319 .329 .762 1.0 .301 .237 .327 .730 .583 1.0 .382 .345 .365 .629 .577 .539 1.0 NOTE: P=7 N=305 ************************************************************************ APPENDIX B: TERMINAL OUTPUT FROM FACAIC OK, SEG FAC_JUL25 ENTER THE NAME OF THE OUTPUT FILE PHY_NO2.OUT THIS PROGRAM WILL COMPUTE THE MAXIMUM LIKELIHOOD PATTERN MATRIX USING INFORMATION CRITERIA (AIC OR CAIC) COPYRIGHT @ JULY 25, 1986 ENTER THE NUMBER OF VARIABLES : <= 24 7 NUMBER OF VARIABLES = 7 ENTER THE NUMBER OF CASES : NO LIMIT 305 NUMBER OF CASES = 305 DO YOU WISH TO USE AIC ENTER Y OR N *****CAPS ONLY***** Y METHOD WILL BE AIC ENTER THE DATA FILE PHY_NO2.DATA THE DATA FILE MAY CONTAIN THE RAW DATA MATRIX OR THE CORRELATION MATRIX IN LOWER DIAGONAL FORM DOES THE DATA FILE CONTAIN THE RAW DATA MATRIX? ENTER Y OR N *****CAPS ONLY***** N DO YOU WANT AN AUTOMATIC SEARCH? ENTER Y OR N *****CAPS ONLY***** Y IS THE DATA IN FREE FORMAT? ENTER Y OR N *****CAPS ONLY***** Y WORKING ON BEST MODEL PROGRAM ENDED NORMALLY **** STOP ************************************************************************ APPENDIX C: FULL OUTPUT FROM FACAIC THIS PROGRAM WILL COMPUTE THE MAXIMUM LIKELIHOOD PATTERN MATRIX USING INFORMATION CRITERIA (AIC OR CAIC) COPYRIGHT @ JULY 25, 1986 NUMBER OF VARIABLES = 7 NUMBER OF CASES = 305 DATA FILE = PHY_NO2.DATA CORRELATION 1 2 3 4 5 6 7 1 1.00000 2 0.80500 1.00000 3 0.85900 0.80100 1.00000 4 0.47300 0.38000 0.43600 1.00000 5 0.39800 0.31900 0.32900 0.76200 1.00000 6 0.30100 0.23700 0.32700 0.73000 0.58300 1.00000 7 0.38200 0.34500 0.36500 0.62900 0.57700 0.53900 1.00000 ************************************************************************ ESTIMATED NUMBER OF FACTORS USING AIC IS 3 ************************************************************************ CHI-SQUARED IS 1.5250 WITH DF = 3 DESCRIPTIVE LEVEL FOR CHI-SQUARE IS 0.6765156 COMMUNALITIES 0.86846 0.75077 0.87528 0.83549 0.75580 0.79234 0.47642 PSI 0.13154 0.24922 0.12471 0.16448 0.24418 0.20765 0.52357 LOADINGS 1 2 3 1 0.89089 0.27344 -0.00245 2 0.84573 0.18846 -0.00023 3 0.90186 0.20417 0.14231 4 0.26391 0.86112 0.15594 5 0.18795 0.84404 -0.08981 6 0.11920 0.71858 0.51164 7 0.24417 0.63772 0.10057 MODEL CORRELATIONS 1 2 3 4 5 6 7 1 0.99999 2 0.80498 0.99999 3 0.85893 0.80117 0.99999 4 0.47019 0.38544 0.43602 0.99997 5 0.39846 0.31804 0.32905 0.76241 0.99998 6 0.30143 0.23611 0.32703 0.73003 0.58296 0.99999 7 0.39166 0.32666 0.36472 0.62927 0.57512 0.53881 0.99999 NORM RESIDUAL CORR 1 2 3 4 5 6 7 1 1.00000 2 0.00011 1.00000 3 0.00053 -0.00096 1.00000 4 0.01907 -0.02688 -0.00012 1.00000 5 -0.00256 0.00389 -0.00031 -0.00207 1.00000 6 -0.00258 0.00389 -0.00017 -0.00014 0.00017 1.00000 7 -0.03679 0.05078 0.00110 -0.00093 0.00526 0.00057 1.00000 RESIDUAL CORRELATION 1 2 3 4 5 6 7 1 -0.00001 2 -0.00002 -0.00001 3 -0.00007 0.00017 -0.00001 4 -0.00281 0.00544 0.00002 -0.00003 5 0.00046 -0.00096 0.00005 0.00041 -0.00002 6 0.00043 -0.00089 0.00003 0.00003 -0.00004 -0.00001 7 0.00966 -0.01834 -0.00028 0.00027 -0.00188 -0.00019 -0.00001 LOG OF DETERMINANT OF MODEL CORR MATRIX IS -5.0437 SCALED LOG-LIKELIHOOD IS 4520.7 NUMBER OF ITERATIONS FOR FINAL MODEL IS 2 MIN TRACE((R^)**(-1)*R)-LOG((R^)**(-1)*R)-P IS 0.51031E-02 NUMBER OF PARAMETERS IS 25.000 PENALTY FOR AIC IS 50.000 AIC IS 4570.7 LOG OF DETERMINANT OF UNRESTRICTED MODEL CORR MATRIX IS -5.0483 SUMMARY OF AIC AND CAIC VALUES: NF NP AIC CAIC DF P-VALUES IER 1 14. 5084.2 5136.3 14. 0.00000 0 2 20. 4580.5 4654.9 8. 0.72759E-02 0 3 25. 4570.7 4663.7 3. 0.67651 0 4 29. 4577.1 4685.0 -1. 39 0 28. 4575.1 4679.3 PROGRAM ENDED NORMALLY ************************************************************************ APPENDIX D: THE FORTRAN PROGRAM FACAIC PROGRAM MAIN 0000001 C DEVELOPED BY HAMPARSUM BOZDOGAN AND DONALD E. RAMIREZ 0000002 C PROGRAMMING ASSISTANTS: LARRY BOBBITT, RENA KAPLAN, AND JAMES SYTA0000003 C MATHEMATICS DEPARTMENT 0000004 C MATHEMATICS-ASTRONOMY BUILDING 0000005 C UNIVERSITY OF VIRGINIA 0000006 C CHARLOTTESVILLE, VIRGINIA 22903 0000007 C COPYRIGHT @ JULY 25, 1986 0000008 INTEGER NROW,NCOL,NCOLSQ,NSYM,MAXVAR 0000009 CHARACTER VERSON*72 0000010 C NOTE: THE MAX NUMBER OF VARIABLES IS MAXVAR 0000011 C THE MAX BLOCK SIZE FOR READING IN DATA IS NROW 0000012 PARAMETER(VERSON='COPYRIGHT @ JULY 25, 1986', 0000013 + MAXVAR=24, 0000014 + NROW=24, 0000015 + NCOL=MAXVAR+1, 0000016 + NCOLSQ=NCOL*NCOL, 0000017 + NSYM=(NCOL*(NCOL+1))/2) 0000018 INTEGER N,P,BESTNF,IOUT,IWORK,DF,NITERS,SETNF,BIGFAC,IEROFC 0000019 DOUBLE PRECISION X,CORR,MEAN,STD,LNLIKE, 0000020 + LAMBDA,PSI,WORK3,FMIN,CHISQ, 0000021 + INVCOR,COMM,RESCOR,OLDPSI, 0000022 + GSCALE,WORK1,WORK2,LNDET,ALPHA, 0000023 + NPARMS,PENALT,AIC,CAIC,CRITER 0000024 LOGICAL DOAIC,RAWDAT,AUTO 0000025 C FOLLOWING IS A LIST OF WORK STORAGE VECTORS AND WHERE THEY 0000026 C ARE USED: 0000027 C WORK1 IS DIMENSIONED AS NCOLSQ=NCOL*NCOL 0000028 C SUBROUTINE FACAIC 0000029 C SUBROUTINE ROTATE 0000030 C SUBROUTINE OFROTA P*(P+1)/2 0000031 C SUBROUTINE PRTOUT 0000032 C SUBROUTINE LINV3F P*P 0000033 C WORK2 IS DIMENSIONED AS NCOLSQ=NCOL*NCOL 0000034 C SUBROUTINE DESRD 0000035 C SUBROUTINE BECOVW 3*P 0000036 C SUBROUTINE FACAIC 0000037 C SUBROUTINE ROTATE 0000038 C SUBROUTINE OFROTA P*P 0000039 C SUBROUTINE PRTOUT P*(P+1)/2 0000040 C WORK3 IS DIMENSIONED AS 5*NCOL 0000041 C SUBROUTINE PRTOUT 0000042 C SUBROUTINE LINV3F P 0000043 C SUBROUTINE FACAIC 0000044 C SUBROUTINE FACTOR 0000045 C SUBROUTINE OFCOMM 5*P 0000046 C SUBROUTINE ROTATE 0000047 C SUBROUTINE OFROTA P 0000048 C IWORK IS DIMENSIONED AS 2*NCOL 0000049 C SUBROUTINE FACAIC 0000050 C SUBROUTINE FACTOR 0000051 C SUBROUTINE OFCOMM 2*P 0000052 DIMENSION X(NROW,NCOL),CORR(NSYM),MEAN(NCOL),STD(NCOL), 0000053 + LAMBDA(NCOL,NCOL),PSI(NCOL), 0000054 + GSCALE(NCOL),WORK2(NCOLSQ),WORK1(NCOLSQ), 0000055 + RESCOR(NSYM),WORK3(NCOL*5),IWORK(2*NCOL), 0000056 + OLDPSI(NCOL,NCOL),CRITER(0:NCOL,5),IEROFC(NCOL), 0000057 + INVCOR(NCOLSQ),COMM(NCOL) 0000058 EXTERNAL START,DESRD,DESOUT,FACAIC,PRTOUT,GETSYM,SYMOUT 0000059 C START OF PROGRAM 0000060 C SETUP INPUT AND OUTPUT FILES 0000061 CALL START(N,P,NROW,NCOL,IOUT,VERSON,DOAIC,RAWDAT,AUTO,SETNF) 0000062 IF (RAWDAT) THEN 0000063 CALL DESRD(N,P,X,MEAN,CORR,STD,IOUT,NROW,WORK2) 0000064 CALL DESOUT(MEAN,CORR,STD,P,IOUT) 0000065 ELSE 0000066 CALL GETSYM(CORR,P,IOUT) 0000067 CALL SYMOUT(CORR,P,IOUT) 0000068 ENDIF 0000069 PRINT 710 0000070 710 FORMAT(/1X,'WORKING ON BEST MODEL') 0000071 CALL FACAIC(N,P,DOAIC,CORR,AUTO,SETNF,BESTNF,ALPHA, 0000072 + NPARMS,PENALT,AIC,CAIC,LNDET,PSI,CRITER,BIGFAC, 0000073 + LAMBDA,INVCOR,COMM,RESCOR,IEROFC,LNLIKE,DF,NITERS, 0000074 + FMIN,CHISQ,OLDPSI,GSCALE,IWORK,WORK1,WORK2,WORK3) 0000075 CALL PRTOUT(N,P,COMM,LNDET,PSI,LAMBDA,BESTNF,RESCOR,IOUT, 0000076 + ALPHA,NPARMS,PENALT,AIC,CAIC,CRITER,BIGFAC,SETNF,AUTO,0000077 + BESTNF,DOAIC,IEROFC,LNLIKE,DF,NITERS,FMIN,CHISQ, 0000078 + CORR,WORK1,WORK2) 0000079 PRINT 720 0000080 WRITE(IOUT,720) 0000081 720 FORMAT(/1X,'PROGRAM ENDED NORMALLY') 0000082 STOP 0000083 END 0000084 ************************************************************************0000085 SUBROUTINE START(N,P,NROW,NCOL,IOUT,VERSON,DOAIC,RAWDAT,AUTO, 0000086 + SETNF) 0000087 INTEGER N,P,NROW,NCOL,IOUT,SETNF 0000088 LOGICAL DOAIC,RAWDAT,AUTO 0000089 CHARACTER VERSON*72 0000090 C LOCAL VARIABLES 0000091 INTEGER IOPT,NIN,UN 0000092 LOGICAL NEW,OLD 0000093 CHARACTER*60 FN 0000094 CHARACTER*80 TITLE 0000095 EXTERNAL QUERY,UGETIO,GTOPFN 0000096 TITLE=' ENTER THE NAME OF THE OUTPUT FILE' 0000097 UN=6 0000098 NEW=.TRUE. 0000099 OLD=.FALSE. 0000100 CALL GTOPFN(UN,FN,NEW,OLD,TITLE) 0000101 IOUT=6 0000102 IOPT=3 0000103 NIN=5 0000104 CALL UGETIO(IOPT,NIN,IOUT) 0000105 IOPT=2 0000106 CALL UGETIO(IOPT,NIN,IOUT) 0000107 PRINT 720, VERSON 0000108 WRITE(IOUT,720) VERSON 0000109 720 FORMAT(/1X,'THIS PROGRAM WILL COMPUTE THE MAXIMUM LIKELIHOOD', 0000110 + ' PATTERN MATRIX',/1X, 0000111 + 'USING INFORMATION CRITERIA (AIC OR CAIC)', 0000112 + /1X,A) 0000113 10 CONTINUE 0000114 PRINT 730,(NCOL-1) 0000115 730 FORMAT(/1X,'ENTER THE NUMBER OF VARIABLES : <= ',I5) 0000116 READ*,P 0000117 IF (P.GT.(NCOL-1)) THEN 0000118 PRINT*,'ERROR IN NUMBER OF VARIABLES' 0000119 GOTO 10 0000120 ENDIF 0000121 PRINT 740,P 0000122 WRITE(IOUT,740) P 0000123 740 FORMAT(/1X,'NUMBER OF VARIABLES = ',I5) 0000124 PRINT 750 0000125 750 FORMAT(/1X,'ENTER THE NUMBER OF CASES : NO LIMIT') 0000126 READ*,N 0000127 PRINT 760,N 0000128 WRITE(IOUT,760) N 0000129 760 FORMAT(/1X,'NUMBER OF CASES = ',I5) 0000130 TITLE='DO YOU WISH TO USE AIC' 0000131 CALL QUERY(TITLE,DOAIC) 0000132 IF (DOAIC) THEN 0000133 PRINT 770 0000134 770 FORMAT(/1X,'METHOD WILL BE AIC') 0000135 ELSE 0000136 PRINT 780 0000137 780 FORMAT(/1X,'METHOD WILL BE CAIC') 0000138 END IF 0000139 UN=5 0000140 NEW=.FALSE. 0000141 OLD=.TRUE. 0000142 TITLE=' ENTER THE DATA FILE' 0000143 CALL GTOPFN(UN,FN,NEW,OLD,TITLE) 0000144 WRITE(IOUT,790) FN 0000145 790 FORMAT(/1X,'DATA FILE = ',A) 0000146 PRINT 800 0000147 800 FORMAT(/1X,'THE DATA FILE MAY CONTAIN THE RAW DATA MATRIX OR THE',0000148 + /1X,'CORRELATION MATRIX IN LOWER DIAGONAL FORM') 0000149 TITLE='DOES THE DATA FILE CONTAIN THE RAW DATA MATRIX?' 0000150 CALL QUERY(TITLE,RAWDAT) 0000151 TITLE='DO YOU WANT AN AUTOMATIC SEARCH?' 0000152 CALL QUERY(TITLE,AUTO) 0000153 IF (.NOT.AUTO) THEN 0000154 PRINT 810 0000155 810 FORMAT(/1X,'ENTER THE NUMBER OF FACTORS') 0000156 READ*,SETNF 0000157 PRINT 820, SETNF 0000158 WRITE(IOUT,820) SETNF 0000159 820 FORMAT(/1X,'NUMBER OF FACTORS HAS BEEN SET EQUAL TO ',I5) 0000160 ENDIF 0000161 RETURN 0000162 END 0000163 ************************************************************************0000164 C SUBROUTINE UGETIO(IOPT,NIN,NOUT) 0000165 C IMSL ROUTINE FOR SETTING AND RETRIEVING INPUT/OUTPUT UNITS 0000166 C END 0000167 ************************************************************************0000168 SUBROUTINE QUERY(TITLE,RESPON) 0000169 LOGICAL RESPON 0000170 CHARACTER*50 TITLE 0000171 C LOCAL VARIABLES 0000172 CHARACTER*1 ANS 0000173 PRINT 710,TITLE 0000174 710 FORMAT(/1X,A) 0000175 10 CONTINUE 0000176 PRINT 720 0000177 720 FORMAT(' ENTER Y OR N *****CAPS ONLY*****') 0000178 READ 730,ANS 0000179 730 FORMAT(A) 0000180 IF((ANS.NE.'Y').AND.(ANS.NE.'N')) THEN 0000181 GOTO 10 0000182 ENDIF 0000183 RESPON=(ANS.EQ.'Y') 0000184 RETURN 0000185 END 0000186 ************************************************************************0000187 SUBROUTINE GTOPFN(UN,FN,NEW,OLD,TITLE) 0000188 C GET AND OPEN FILENAME 0000189 C REVISED VERSION FEB 2 1985 0000190 INTEGER UN 0000191 LOGICAL NEW,OLD 0000192 CHARACTER*60 FN 0000193 CHARACTER*80 TITLE 0000194 C UN (INPUT) UNIT TO OPEN FILE ON 0000195 C FN (OUTPUT) NAME OF FILE THAT HAS BEEN OPENED 0000196 C NEW (INPUT) 0000197 C .TRUE. IF THE FILE IS REQUIRED TO BE A NEW FILE 0000198 C .FALSE. OTHERWISE 0000199 C OLD (INPUT) 0000200 C .TRUE. IF THE FILE IS REQUIRED TO BE A OLD FILE 0000201 C .FALSE. OTHERWISE 0000202 C TITLE (INPUT) PROMPT TO BE PRINTED WHEN ASKING FOR FN 0000203 C LOCAL VARIABLES 0000204 LOGICAL EX 0000205 10 CONTINUE 0000206 PRINT* 0000207 PRINT 710,TITLE 0000208 710 FORMAT(A) 0000209 READ 710,FN 0000210 INQUIRE(FILE=FN,EXIST=EX) 0000211 IF(NEW.AND.EX) THEN 0000212 PRINT*,'FILE = ' 0000213 PRINT*,' ',FN 0000214 PRINT*,'EXISTS. PLEASE TRY AGAIN.' 0000215 GOTO 10 0000216 ENDIF 0000217 IF(OLD.AND.(.NOT.EX)) THEN 0000218 PRINT*,' FILE = ' 0000219 PRINT*,' ',FN 0000220 PRINT*,' DOES NOT EXIST. PLEASE TRY AGAIN.' 0000221 GOTO 10 0000222 ENDIF 0000223 OPEN(UNIT=UN,FILE=FN,ERR=20) 0000224 REWIND UN 0000225 RETURN 0000226 20 CONTINUE 0000227 PRINT*,' ERROR IN OPENING FILE =' 0000228 PRINT*,FN 0000229 PRINT*,' PLEASE TRY AGAIN.' 0000230 GOTO 10 0000231 END 0000232 ************************************************************************0000233 SUBROUTINE GETSYM(X,P,IOUT) 0000234 INTEGER P,IOUT 0000235 DOUBLE PRECISION X(P*(P+1)/2) 0000236 C LOCAL VARIABLES 0000237 INTEGER I,INUNIT 0000238 LOGICAL FIRST,FREE 0000239 CHARACTER FMT*80,TITLE*50 0000240 EXTERNAL QUERY 0000241 SAVE FIRST,FMT 0000242 DATA FIRST/.TRUE./ 0000243 INUNIT=5 0000244 TITLE='IS THE DATA IN FREE FORMAT?' 0000245 CALL QUERY(TITLE,FREE) 0000246 IF (FIRST.AND.(.NOT.FREE)) THEN 0000247 PRINT 710 0000248 710 FORMAT(/'ENTER THE FORTRAN FORMAT FOR THE DATA') 0000249 READ 720, FMT 0000250 720 FORMAT(A) 0000251 FIRST=.FALSE. 0000252 ENDIF 0000253 IF(.NOT.FREE) THEN 0000254 READ(INUNIT,FMT) (X(I),I=1,P*(P+1)/2) 0000255 ELSE 0000256 READ(INUNIT,*,ERR=30) (X(I),I=1,P*(P+1)/2) 0000257 ENDIF 0000258 RETURN 0000259 30 CONTINUE 0000260 PRINT 730,I 0000261 WRITE(IOUT,730) I 0000262 730 FORMAT(/1X,'ERROR IN READING DATA MATRIX' 0000263 + /1X,'ERROR OCCURRED WHILE READING CASE ',I5) 0000264 STOP 0000265 END 0000266 ************************************************************************0000267 SUBROUTINE SYMOUT(X,P,IOUT) 0000268 INTEGER P,IOUT 0000269 DOUBLE PRECISION X(P*(P+1)/2) 0000270 C LOCAL VARIABLES 0000271 INTEGER IOPT,NC 0000272 CHARACTER*20 ITITLE 0000273 EXTERNAL USWSM 0000274 WRITE(IOUT,710) 0000275 710 FORMAT(/) 0000276 ITITLE='CORRELATION' 0000277 NC=20 0000278 IOPT=2 0000279 CALL USWSM(ITITLE,NC,X,P,IOPT) 0000280 WRITE(IOUT,710) 0000281 RETURN 0000282 END 0000283 ************************************************************************0000284 SUBROUTINE DESRD(N,P,X,MEAN,CORR,STD,IOUT,NROW,WORK2) 0000285 INTEGER N,P,NROW,IOUT 0000286 DOUBLE PRECISION X(NROW,P),STD(P), 0000287 + MEAN(P),CORR(P*(P+1)/2),WORK2(3*P) 0000288 C LOCAL VARIABLES 0000289 DOUBLE PRECISION XMISS,ONE,WT 0000290 INTEGER I,J,NLEFT,IMAT,INCD,IER,IN,INUNIT,INDSYM,NBR(6) 0000291 LOGICAL FREE 0000292 CHARACTER TITLE*50 0000293 PARAMETER(ONE=1.0D0) 0000294 INTRINSIC REAL,NINT,MIN 0000295 EXTERNAL QUERY,GETDAT,BECOVW,INDSYM 0000296 TITLE='IS THE DATA IN FREE FORMAT?' 0000297 CALL QUERY(TITLE,FREE) 0000298 NBR(1)=P 0000299 NBR(2)=N 0000300 NBR(3)=NINT(MIN(REAL(NROW),REAL(N) )) 0000301 NBR(5)=0 0000302 NBR(6)=-1 0000303 WT=-ONE 0000304 XMISS=0.0 0000305 INCD=0 0000306 IMAT=0 0000307 NLEFT=N 0000308 INUNIT=5 0000309 10 CONTINUE 0000310 IMAT=IMAT+1 0000311 IN=NINT( MIN( REAL(NLEFT) , REAL(NROW) ) ) 0000312 CALL GETDAT(FREE,X,NROW,P,IN,INUNIT,IOUT) 0000313 IF(IMAT.EQ.1) THEN 0000314 WRITE(IOUT,710) 'FIRST CASE IS:' 0000315 WRITE(IOUT,720) (X(1,J),J=1,P) 0000316 ENDIF 0000317 NBR(4)=IMAT 0000318 CALL BECOVW(X,NROW,WT,NBR,XMISS,MEAN,CORR,INCD,WORK2,IER) 0000319 NLEFT=NLEFT-IN 0000320 IF(NLEFT.GT.0) GOTO 10 0000321 REWIND 5 0000322 WRITE(IOUT,710) 'LAST CASE IS:' 0000323 710 FORMAT(/1X,A) 0000324 WRITE(IOUT,720) (X(IN,J),J=1,P) 0000325 720 FORMAT(1X,6G12.4) 0000326 DO 20,I=1,P 0000327 J=INDSYM(I,I) 0000328 STD(I)=CORR(J) 0000329 CORR(J)=ONE 0000330 20 CONTINUE 0000331 RETURN 0000332 END 0000333 ************************************************************************0000334 C SUBROUTINE BECOVW(X,IX,WT,NBR,XMISS,XM,VCV,INCD,WK,IER) 0000335 C IMSL ROUTINE FOR COMPUTING CORRELATION IN PIECES 0000336 C TERMINAL ERRORS: 0000337 C IER=129 INDICATES THAT NBR(3)*(NBR(4)-1) EXCEEDS NBR(2) 0000338 C (WHEN NBR(4) IS POSITIVE). 0000339 C IER=130 INDICATES THAT NBR(1) IS LESS THAN 1 OR NBR(2) IS LESS 0000340 C THAN 2 OR NBR(3) EXCEEDS NBR(2). 0000341 C IER=131 INDICATES THAT NBR(4)=1 AND WT(1) IS NONNEGATIVE 0000342 C (IMPLYING WEIGHTS ARE TO BE USED) AND THAT SOME WEIGHTS0000343 C ARE NEGATIVE. 0000344 C END 0000345 ************************************************************************0000346 SUBROUTINE GETDAT(FREE,X,NROW,P,IN,INUNIT,IOUT) 0000347 INTEGER NROW,P,IN,INUNIT,IOUT 0000348 DOUBLE PRECISION X(NROW,P) 0000349 LOGICAL FREE 0000350 C LOCAL VARIABLES 0000351 INTEGER I,J 0000352 LOGICAL FIRST 0000353 CHARACTER*80 FMT 0000354 SAVE FIRST,FMT 0000355 DATA FIRST/.TRUE./ 0000356 IF(FIRST.AND.(.NOT.FREE)) THEN 0000357 PRINT 710 0000358 710 FORMAT(/'ENTER THE FORTRAN FORMAT FOR THE DATA') 0000359 READ 720, FMT 0000360 720 FORMAT(A) 0000361 FIRST=.FALSE. 0000362 ENDIF 0000363 IF(.NOT.FREE) THEN 0000364 DO 10,I=1,IN 0000365 READ(INUNIT,FMT) (X(I,J),J=1,P) 0000366 10 CONTINUE 0000367 ELSE 0000368 DO 20,I=1,IN 0000369 READ(INUNIT,*,ERR=30) (X(I,J),J=1,P) 0000370 20 CONTINUE 0000371 ENDIF 0000372 RETURN 0000373 30 CONTINUE 0000374 PRINT 730,I 0000375 WRITE(IOUT,730) I 0000376 730 FORMAT(/1X,'ERROR IN READING DATA MATRIX' 0000377 + /1X,'ERROR OCCURRED WHILE READING CASE ',I5) 0000378 STOP 0000379 END 0000380 ************************************************************************0000381 FUNCTION INDSYM(I,J) 0000382 INTEGER INDSYM,I,J 0000383 INDSYM=I*(I-1)/2 + J 0000384 RETURN 0000385 END 0000386 ************************************************************************0000387 SUBROUTINE DESOUT(MEAN,CORR,STD,P,IOUT) 0000388 INTEGER P,IOUT 0000389 DOUBLE PRECISION MEAN(P),CORR(P*(P+1)/2),STD(P) 0000390 C LOCAL VARIABLES 0000391 INTEGER INC,IOPT,NC 0000392 CHARACTER*20 ITITLE 0000393 EXTERNAL USWFV,USWSM 0000394 WRITE(IOUT,710) 0000395 710 FORMAT(/) 0000396 ITITLE='MEANS OF DATA MATRIX' 0000397 NC=20 0000398 INC=1 0000399 IOPT=2 0000400 CALL USWFV(ITITLE,NC,MEAN,P,INC,IOPT) 0000401 WRITE(IOUT,710) 0000402 ITITLE='STANDARD DEVIATIONS' 0000403 NC=20 0000404 INC=1 0000405 IOPT=2 0000406 CALL USWFV(ITITLE,NC,STD,P,INC,IOPT) 0000407 WRITE(IOUT,710) 0000408 ITITLE='CORRELATION' 0000409 NC=20 0000410 IOPT=2 0000411 CALL USWSM(ITITLE,NC,CORR,P,IOPT) 0000412 WRITE(IOUT,710) 0000413 RETURN 0000414 END 0000415 ************************************************************************0000416 C SUBROUTINE USWFV(ITITLE,NC,A,M,INC,IOPT) 0000417 C IMSL ROUTINE FOR PRINTING A VECTOR 0000418 C END 0000419 ************************************************************************0000420 C SUBROUTINE USWSM(ITITLE,NC,A,M,IOPT) 0000421 C IMSL ROUTINE FOR PRINTING A SYMMETRIC MATRIX 0000422 C END 0000423 ************************************************************************0000424 SUBROUTINE PRTOUT(N,P,COMM,LNDET,PSI,LAMBDA,NF,RESCOR,IOUT,ALPHA, 0000425 + NPARMS,PENALT,AIC,CAIC,CRITER,BIGFAC,SETNF,AUTO,0000426 + BESTNF,DOAIC,IEROFC,LNLIKE,DF,NITERS,FMIN,CHISQ,0000427 + CORR,WORK1,WORK2) 0000428 INTEGER N,P,NF,IOUT,IEROFC(P),BESTNF,DF,NITERS,SETNF,BIGFAC 0000429 DOUBLE PRECISION PSI(P),LAMBDA(P,P),ALPHA,LNDET,COMM(P),LNLIKE, 0000430 + FMIN,CHISQ,RESCOR((P*(P+1))/2),WORK1(P,P), 0000431 + NPARMS,PENALT,AIC,CAIC,CRITER(0:P,5), 0000432 + CORR((P*(P+1))/2),WORK2((P*(P+1))/2) 0000433 LOGICAL DOAIC,AUTO 0000434 C LOCAL VARIABLES 0000435 INTEGER NC,INC,IOPT,IJOB,I,LASTNF,IER 0000436 DOUBLE PRECISION B,D1,D2,DET 0000437 CHARACTER ITITLE*20 0000438 INTRINSIC LOG 0000439 EXTERNAL USWFV,USWFM,VMULFP,LINV3F,USWSM,VCVTFS 0000440 IF (DOAIC) THEN 0000441 ITITLE='AIC' 0000442 ELSE 0000443 ITITLE='CAIC' 0000444 END IF 0000445 WRITE(IOUT,720) ITITLE,NF 0000446 720 FORMAT(/1X,72('*'),//1X,'ESTIMATED NUMBER OF FACTORS USING ' 0000447 + ,A,' IS ',I5,//1X,72('*'),/) 0000448 IF (BESTNF.GT.0) THEN 0000449 WRITE(IOUT,730) CHISQ,DF 0000450 730 FORMAT(/1X,'CHI-SQUARED IS ',G16.5,' WITH DF = ',I5) 0000451 IF (DF.GT.0) THEN 0000452 WRITE(IOUT,740) ALPHA 0000453 740 FORMAT(/1X,'DESCRIPTIVE LEVEL FOR CHI-SQUARE IS ',F16.7) 0000454 ELSE 0000455 WRITE(IOUT,750) 0000456 750 FORMAT(/1X,'DESCRIPTIVE LEVEL FOR CHI-SQUARE', 0000457 + ' CANNOT BE COMPUTED') 0000458 END IF 0000459 ITITLE = 'COMMUNALITIES' 0000460 NC = 20 0000461 INC = 1 0000462 IOPT = 2 0000463 WRITE(IOUT,760) 0000464 760 FORMAT(1X) 0000465 CALL USWFV(ITITLE,NC,COMM,P,INC,IOPT) 0000466 WRITE(IOUT,760) 0000467 C760 FORMAT(1X) 0000468 ITITLE='PSI' 0000469 NC=20 0000470 INC=1 0000471 IOPT=2 0000472 CALL USWFV(ITITLE,NC,PSI,P,INC,IOPT) 0000473 WRITE(IOUT,760) 0000474 C760 FORMAT(1X) 0000475 ITITLE='LOADINGS' 0000476 NC=20 0000477 IOPT=2 0000478 CALL USWFM(ITITLE,NC,LAMBDA,P,P,NF,IOPT) 0000479 C NOTE: WORK1 IS THE CORRELATION MATRIX UNDER 0000480 C THE CURRENT FACTOR MODEL 0000481 CALL VMULFP(LAMBDA,LAMBDA,P,NF,P,P,P,WORK1,P,IER) 0000482 DO 10,I=1,P 0000483 WORK1(I,I)=WORK1(I,I)+PSI(I) 0000484 10 CONTINUE 0000485 C PUT MODEL CORRELATION MATRIX INTO SYMMETRIC STORAGE IN WORK2 0000486 CALL VCVTFS(WORK1,P,P,WORK2) 0000487 ITITLE='MODEL CORRELATIONS' 0000488 NC=20 0000489 IOPT=2 0000490 WRITE(IOUT,760) 0000491 C760 FORMAT(1X) 0000492 CALL USWSM(ITITLE,NC,WORK2,P,IOPT) 0000493 WRITE(IOUT,760) 0000494 C760 FORMAT(1X) 0000495 ITITLE='NORM RESIDUAL CORR' 0000496 NC=20 0000497 IOPT=2 0000498 CALL USWSM(ITITLE,NC,RESCOR,P,IOPT) 0000499 C COMPUTE RESIDUAL CORRELATION MATRIX 0000500 DO 20,I=1,(P*(P+1))/2 0000501 WORK2(I)=WORK2(I)-CORR(I) 0000502 20 CONTINUE 0000503 WRITE(IOUT,760) 0000504 C760 FORMAT(1X) 0000505 ITITLE='RESIDUAL CORRELATION' 0000506 NC=20 0000507 IOPT=2 0000508 CALL USWSM(ITITLE,NC,WORK2,P,IOPT) 0000509 C COMPUTE THE DETERMINANT OF THE MODEL CORRELATION MATRIX WORK1 0000510 IJOB=4 0000511 B=0.0D0 0000512 D1=0.0D0 0000513 CALL LINV3F(WORK1,B,IJOB,P,P,D1,D2,WORK2,IER) 0000514 DET=D1*2.0D0**D2 0000515 WRITE(IOUT,770) LOG(DET) 0000516 770 FORMAT(/1X,'LOG OF DETERMINANT OF MODEL CORR MATRIX IS ', 0000517 + G16.5) 0000518 WRITE(IOUT,780) LNLIKE 0000519 780 FORMAT(/1X,'SCALED LOG-LIKELIHOOD IS ',G16.5) 0000520 WRITE(IOUT,790) NITERS 0000521 790 FORMAT(/1X,'NUMBER OF ITERATIONS FOR FINAL MODEL IS ',I5) 0000522 WRITE(IOUT,800) FMIN 0000523 800 FORMAT(/1X,'MIN TRACE((R^)**(-1)*R)-LOG((R^)**(-1)*R)-P IS ', 0000524 + G16.5) 0000525 WRITE(IOUT,810) NPARMS 0000526 810 FORMAT(/1X,'NUMBER OF PARAMETERS IS ',G16.5) 0000527 IF (DOAIC) THEN 0000528 WRITE(IOUT,820) PENALT 0000529 820 FORMAT(/1X,'PENALTY FOR AIC IS ',G16.5) 0000530 WRITE(IOUT,830) AIC 0000531 830 FORMAT(/1X,'AIC IS ',G16.5) 0000532 ELSE 0000533 WRITE(IOUT,820) PENALT 0000534 C820 FORMAT(/1X,'PENALTY FOR CAIC IS ',G16.5) 0000535 WRITE(IOUT,830) CAIC 0000536 C830 FORMAT(/1X,'CAIC IS ',G16.5) 0000537 ENDIF 0000538 ENDIF 0000539 WRITE(IOUT,840) LNDET 0000540 840 FORMAT(/1X,'LOG OF DETERMINANT OF UNRESTRICTED MODEL ' 0000541 + ,'CORR MATRIX IS ',G16.5) 0000542 IF (AUTO) THEN 0000543 LASTNF=BIGFAC 0000544 ELSE 0000545 LASTNF=SETNF 0000546 ENDIF 0000547 WRITE(IOUT,850) 0000548 850 FORMAT(/1X,'SUMMARY OF AIC AND CAIC VALUES:') 0000549 WRITE(IOUT,860) 'NF','NP','AIC','CAIC','DF','P-VALUES','IER' 0000550 860 FORMAT(/1X,A4,1X,A4,1X,A10,7X,A10,7X,A4,1X,A16,1X,A4) 0000551 DO 30,I=1,LASTNF 0000552 IF (CRITER(I,4).GT.0.5D0) THEN 0000553 WRITE(IOUT,870) I,CRITER(I,1),CRITER(I,2),CRITER(I,3), 0000554 + CRITER(I,4),CRITER(I,5),IEROFC(I) 0000555 870 FORMAT(/1X,I4,1X,F4.0,1X,G16.5,1X,G16.5,1X,F4.0,1X,G16.5, 0000556 + 1X,I4) 0000557 ELSE 0000558 WRITE(IOUT,880) I,CRITER(I,1),CRITER(I,2),CRITER(I,3), 0000559 + CRITER(I,4),IEROFC(I) 0000560 880 FORMAT(/1X,I4,1X,F4.0,1X,G16.5,1X,G16.5,1X,F4.0,1X,16X, 0000561 + 1X,I4) 0000562 ENDIF 0000563 30 CONTINUE 0000564 WRITE(IOUT,890) 0,CRITER(0,1),CRITER(0,2),CRITER(0,3) 0000565 890 FORMAT(/1X,I4,1X,F4.0,1X,G16.5,1X,G16.5) 0000566 RETURN 0000567 END 0000568 ************************************************************************0000569 C SUBROUTINE LINV3F(A,B,IJOB,N,N,IA,D1,D2,WKAREA,IER) 0000570 C IMSL ROUTINE FOR COMPUTING DETERMINANTS 0000571 C TERMINAL ERROR: 0000572 C IER=130 INDICATES THAT MATRIX A IS ALGORITHMICALLY SINGULAR. 0000573 C END 0000574 ************************************************************************0000575 C SUBROUTINE VCVTFS(A,N,IA,B) 0000576 C IMSL ROUTINE TO CONVERT A MATRIX FROM FULL TO SYMMETRIC STORAGE 0000577 C END 0000578 ************************************************************************0000579 C SUBROUTINE VMULFP(A,B,L,M,N,IA,IB,C,IC,IER) 0000580 C IMSL ROUTINE FOR MULTIPLICATION OF A BY B-TRANSPOSE 0000581 C TERMINAL ERROR: 0000582 C IER=129 INDICATES A, B, OR C WAS DIMENSIONED INCORRECTLY. 0000583 C END 0000584 ************************************************************************0000585 SUBROUTINE FACAIC(N,P,DOAIC,CORR,AUTO,SETNF,BESTNF,ALPHA, 0000586 + NPARMS,PENALT,AIC,CAIC,LNDET,PSI,CRITER,BIGFAC, 0000587 + LAMBDA,INVCOR,COMM,RESCOR,IEROFC,LNLIKE,DF, 0000588 + NITERS,FMIN,CHISQ,OLDPSI,GSCALE, 0000589 + IWORK,WORK1,WORK2,WORK3) 0000590 C DEVELOPED BY HAMPARSUM BOZDOGAN AND DONALD E. RAMIREZ 0000591 C PROGRAMMING ASSISTANTS: LARRY BOBBITT, RENA KAPLAN, AND JAMES SYTA0000592 C MATHEMATICS DEPARTMENT 0000593 C MATHEMATICS-ASTRONOMY BUILDING 0000594 C UNIVERSITY OF VIRGINIA 0000595 C CHARLOTTESVILLE, VIRGINIA 22903 0000596 C COPYRIGHT @ JULY 25, 1986 0000597 C N INPUT : INTEGER 0000598 C NUMBER OF CASES 0000599 C P INPUT : INTEGER 0000600 C NUMBER OF VARIABLES 0000601 C DOAIC INPUT : LOGICAL 0000602 C .TRUE. MEANS AIC CRITERION 0000603 C .FALSE. MEANS CAIC CRITERION 0000604 C CORR INPUT : DOUBLE PRECISION 0000605 C DIMENSION (P*(P+1)/2) 0000606 C CORRELATION MATRIX IN SYMMETRIC STORAGE 0000607 C AUTO INPUT : LOGICAL 0000608 C .TRUE. INDICATES AN AUTOMATIC SEARCH FOR THE BEST 0000609 C NUMBER OF FACTORS 0000610 C .FALSE. INDICATES THAT THE USER HAS SUPPLIED THE 0000611 C NUMBER OF FACTORS IN SETNF 0000612 C SETNF INPUT : INTEGER 0000613 C USER SPECIFIED NUMBER OF FACTORS WHEN AUTO=.FALSE. 0000614 C BESTNF OUTPUT : INTEGER 0000615 C BEST NUMBER OF FACTORS 0000616 C ALPHA OUTPUT : DOUBLE PRECISION 0000617 C P-VALUE FOR GOODNESS-OF-FIT USING BESTNF FACTORS 0000618 C NPARMS OUTPUT : DOUBLE PRECISION 0000619 C NUMBER OF PARAMETERS 0000620 C PENALT OUTPUT : DOUBLE PRECISION 0000621 C PENALTY COMPONENT FOR AIC OR CAIC 0000622 C AIC OUTPUT : DOUBLE PRECISION 0000623 C AKAIKE'S INFORMATION CRITERION FOR MODEL SELECTION 0000624 C CAIC OUTPUT : DOUBLE PRECISION 0000625 C CONSISTENT AIC FOR MODEL SELECTION 0000626 C LNDET OUTPUT : DOUBLE PRECISION 0000627 C LOG OF DETERMINANT OF CORRELATION MATRIX 0000628 C PSI OUTPUT : DOUBLE PRECISION 0000629 C DIMENSION (P) 0000630 C PSI VECTOR OF ERROR VARIANCES 0000631 C CRITER OUTPUT: DOUBLE PRECISION 0000632 C DIMENSION (0:P,5) 0000633 C MATRIX OF NUMBER OF PARAMETERS, AIC VALUES, 0000634 C CAIC VALUES, DEGREES-OF-FREEDOM, AND P-VALUES 0000635 C LAMBDA OUTPUT : DOUBLE PRECISION 0000636 C DIMENSION (P,P) 0000637 C LAMBDA MATRIX OF FACTOR LOADINGS 0000638 C INVCOR OUTPUT : DOUBLE PRECISION 0000639 C DIMENSION (P*(P+1)/2) 0000640 C INVERSE OF CORRELATION MATRIX IN SYMMETRIC STORAGE0000641 C COMM OUTPUT : DOUBLE PRECISION 0000642 C DIMENSION (P) 0000643 C VECTOR OF COMMUNALITIES 0000644 C RESCOR OUTPUT : DOUBLE PRECISION 0000645 C DIMENSION (P*(P+1)/2) 0000646 C NORMALIZED RESIDUAL CORRELATION MATRIX IN 0000647 C SYMMETRIC STORAGE 0000648 C (I,J) OFF-DIAGONAL ELEMENT IS EQUAL TO 0000649 C (SAMPLE CORR(I,J)-MODEL CORR(I,J))/SQRT(PSI(I)*PSI(J))0000650 C IEROFC OUTPUT : INTEGER 0000651 C DIMENSION (0:P) 0000652 C IMSL ERROR FLAG FOR OFCOMM FOR EACH MODEL 0000653 C LNLIKE OUTPUT : DOUBLE PRECISION 0000654 C SCALED LOG-LIKELIHOOD 0000655 C DF OUTPUT : INTEGER 0000656 C DEGREES OF FREEDOM IN FINAL MODEL 0000657 C NITERS OUTPUT : INTEGER 0000658 C NUMBER OF ITERATIONS USED IN OFCOMM IN FINAL MODEL0000659 C FMIN OUTPUT : DOUBLE PRECISION 0000660 C VALUE OF THE FUNCTION MINIMUM UNDER THE MAXIMUM 0000661 C LIKELIHOOD PROCEDURE IN OFCOMM IN FINAL MODEL 0000662 C CHISQ OUTPUT : DOUBLE PRECISION 0000663 C THE CHI-SQUARED STATISTIC IN OFCOMM IN FINAL MODEL 0000664 C OLDPSI WORK : DOUBLE PRECISION 0000665 C DIMENSION (P,P) 0000666 C STORES PREVIOUS PSI VALUES FOR NEXT CALL TO OFCOMM 0000667 C GSCALE WORK : DOUBLE PRECISION 0000668 C DIMENSION (P) 0000669 C SCALING VECTOR NEEDED BY SOME IMSL ROTATION PROGRAMS0000670 C IWORK WORK : INTEGER 0000671 C DIMENSION (2*P) 0000672 C WORK VECTOR USED BY OFCOMM 0000673 C WORK1 WORK : DOUBLE PRECISION 0000674 C DIMENSION (P*(P+1)/2) 0000675 C WORK VECTOR USED BY OFROTA AND LUDECP 0000676 C WORK2 WORK : DOUBLE PRECISION 0000677 C DIMENSION (P,P) 0000678 C WORK MATRIX USED BY OFROTA 0000679 C WORK3 WORK : DOUBLE PRECISION 0000680 C DIMENSION (5*P) 0000681 C WORK VECTOR USED BY OFCOMM AND OFROTA 0000682 INTEGER N,P,BESTNF,IWORK(P*2),IEROFC(0:P),DF,NITERS,SETNF,BIGFAC 0000683 DOUBLE PRECISION LNLIKE,LAMBDA(P,P),PSI(P),FMIN,CHISQ, 0000684 + NPARMS,PENALT,AIC,CAIC,CRITER(0:P,5), 0000685 + INVCOR(P*(P+1)/2),COMM(P),RESCOR(P*(P+1)/2), 0000686 + OLDPSI(P,P),ALPHA,GSCALE(P), 0000687 + WORK1(P*(P+1)/2),WORK2(P,P),WORK3(P*5), 0000688 + CORR(P*(P+1)/2),LNDET 0000689 LOGICAL DOAIC,AUTO 0000690 C LOCAL VARIABLES 0000691 INTEGER NF,IZERO,LEVEL,LEVOLD,IER 0000692 DOUBLE PRECISION D1,D2,LN2PI,DET,ZERLIK,XBIGFC 0000693 LOGICAL FIRST 0000694 PARAMETER (IZERO=0) 0000695 INTRINSIC SQRT,INT,LOG,ATAN,MIN,DBLE 0000696 EXTERNAL FACTOR,DOCRIT,ROTATE,LUDECP,UERSET 0000697 C FOR LEVEL=1 IMSL MESSAGES ARE PRINTED ONLY IF IER>128 0000698 LEVEL=1 0000699 C SET IMSL MESSAGE LEVEL 0000700 CALL UERSET(LEVEL,LEVOLD) 0000701 XBIGFC=(2.0D0*P+1.0D0-SQRT(8.0D0*P+1.0D0))/2.0D0 0000702 C NEED LARGEST INTEGER GREATER THAN OF EQUAL TO XBIGFC 0000703 BIGFAC=INT(XBIGFC+.999999) 0000704 FIRST=.TRUE. 0000705 LN2PI=LOG(8.0D0*ATAN(1.0D0)) 0000706 C COMPUTING DETERMINANT OF CORR 0000707 CALL LUDECP(CORR,WORK1,P,D1,D2,IER) 0000708 DET=D1*2.0D0**D2 0000709 LNDET=LOG(DET) 0000710 C ZERLIK IS -2*LOG(LIKELIHOOD) IN THE UNRESTRICTED MODEL 0000711 ZERLIK=N*(P+LNDET+LN2PI*P) 0000712 C COMPUTING INFORMATION CRITERION 0000713 CALL DOCRIT(N,P,IZERO,ZERLIK,BESTNF,DOAIC,NPARMS,PENALT,AIC,CAIC, 0000714 + CRITER) 0000715 IF (AUTO) THEN 0000716 DO 20,NF=1,BIGFAC 0000717 C COMPUTE FACTOR LOADINGS 0000718 CALL FACTOR(N,P,NF,CORR,LNDET,PSI,LAMBDA,INVCOR,COMM, 0000719 + RESCOR,LNLIKE,DF,NITERS,FMIN,CHISQ,FIRST,ALPHA, 0000720 + IEROFC,GSCALE,OLDPSI,IWORK,WORK3) 0000721 C IWORK(1) IS DEGREES-OF-FREEDOM 0000722 CRITER(NF,4)=DBLE(IWORK(1)) 0000723 CRITER(NF,5)=ALPHA 0000724 C COMPUTING INFORMATION CRITERION 0000725 CALL DOCRIT(N,P,NF,LNLIKE,BESTNF,DOAIC,NPARMS,PENALT,AIC, 0000726 + CAIC,CRITER) 0000727 20 CONTINUE 0000728 ELSE 0000729 IF (SETNF.EQ.IZERO) RETURN 0000730 SETNF=MIN(SETNF,BIGFAC) 0000731 DO 30,NF=1,SETNF 0000732 C COMPUTE FACTOR LOADINGS 0000733 CALL FACTOR(N,P,NF,CORR,LNDET,PSI,LAMBDA,INVCOR,COMM, 0000734 + RESCOR,LNLIKE,DF,NITERS,FMIN,CHISQ,FIRST,ALPHA, 0000735 + IEROFC,GSCALE,OLDPSI,IWORK,WORK3) 0000736 C IWORK(1) IS DEGREES-OF-FREEDOM 0000737 CRITER(NF,4)=DBLE(IWORK(1)) 0000738 CRITER(NF,5)=ALPHA 0000739 C COMPUTING INFORMATION CRITERION 0000740 CALL DOCRIT(N,P,NF,LNLIKE,BESTNF,DOAIC,NPARMS,PENALT,AIC, 0000741 + CAIC,CRITER) 0000742 30 CONTINUE 0000743 BESTNF=SETNF 0000744 ENDIF 0000745 IF (BESTNF.GT.IZERO) THEN 0000746 FIRST=.FALSE. 0000747 C COMPUTING FACTOR MODEL WITH NF=BESTNF 0000748 CALL FACTOR(N,P,BESTNF,CORR,LNDET,PSI,LAMBDA,INVCOR,COMM, 0000749 + RESCOR,LNLIKE,DF,NITERS,FMIN,CHISQ,FIRST,ALPHA, 0000750 + IEROFC,GSCALE,OLDPSI,IWORK,WORK3) 0000751 C COMPUTE INFORMATION CRITERION 0000752 CALL DOCRIT(N,P,BESTNF,LNLIKE,BESTNF,DOAIC,NPARMS,PENALT,AIC, 0000753 + CAIC,CRITER) 0000754 C COMPUTE ROTATED FACTOR LOADINGS 0000755 CALL ROTATE(P,BESTNF,LAMBDA,WORK1,WORK2,WORK3) 0000756 ENDIF 0000757 RETURN 0000758 END 0000759 ************************************************************************0000760 C SUBROUTINE LUDECP(A,UL,N,D1,D2,IER) 0000761 C IMSL ROUTINE FOR COMPUTING DETERMINANTS OF SYMMETRIC MATRICES 0000762 C TERMINAL ERROR: 0000763 C IER=129 INDICATES THAT MATRIX A IS ALGORITHMICALLY NOT POSITIVE0000764 C DEFINITE. 0000765 C END 0000766 ************************************************************************0000767 C SUBROUTINE UERSET(LEVEL,LEVOLD) 0000768 C IMSL ROUTINE FOR SUPPRESSION OF WARNING ERRORS 0000769 C END 0000770 ************************************************************************0000771 SUBROUTINE FACTOR(N,P,NF,CORR,LNDET,PSI,LAMBDA,INVCOR,COMM, 0000772 + RESCOR,LNLIKE,DF,NITERS,FMIN,CHISQ,FIRST,ALPHA, 0000773 + IEROFC,GSCALE,OLDPSI,IWORK,WORK3) 0000774 INTEGER N,P,NF,IWORK(P*2),IEROFC(P),DF,NITERS 0000775 DOUBLE PRECISION LNLIKE,GSCALE(P),RESCOR(P*(P+1)/2),COMM(P), 0000776 + INVCOR(P*(P+1)/2),OLDPSI(P,P),LNDET,ALPHA, 0000777 + CORR(P*(P+1)/2),PSI(P),LAMBDA(P,P),WORK3(P*5), 0000778 + FMIN,CHISQ 0000779 LOGICAL FIRST 0000780 C LOCAL VARIABLES 0000781 INTEGER IND,IV,I,MAXIT,MAXTRY 0000782 DOUBLE PRECISION EPS,EPSE,LN2PI 0000783 INTRINSIC DBLE,ATAN,LOG 0000784 EXTERNAL OFCOMM 0000785 C CONVERGENCE CRITERIA USED BY THE IMSL ROUTINE OFCOMM: 0000786 C MAXIT: MAXIMUM NUMBER OF ITERATIONS 0000787 C IMSL SUGGESTS 30 0000788 C MAXTRY: MAXIMUM NUMBER OF SEARCHES FOR A DESCENT POINT 0000789 C IMSL SUGGESTS 25 0000790 C EPS: STOPPING CRITERION 0000791 C IMSL SUGGESTS 0.005 0000792 C EPSE: APPROXIMATE SECOND DERIVATIVE CRITERION 0000793 C IMSL SUGGESTS 0.1 0000794 PARAMETER(MAXIT=120, 0000795 + MAXTRY=100, 0000796 + EPS=0.0050D0, 0000797 + EPSE=0.1D0) 0000798 DATA IV/0/ 0000799 LN2PI=LOG(8.0D0*ATAN(1.0D0)) 0000800 IF (.NOT.FIRST) THEN 0000801 DO 10,I=1,P 0000802 PSI(I)=OLDPSI(I,NF) 0000803 10 CONTINUE 0000804 ENDIF 0000805 C IND=3 IMPLIES MAX-LIKELIHOOD PROCEDURE 0000806 IND=3 0000807 ALPHA=-99.0D0 0000808 C COMPUTING UNROTATED FACTOR LOADING MATRIX 0000809 CALL OFCOMM(CORR,P,NF,IND,N,IV,MAXIT,MAXTRY,EPS,EPSE, 0000810 + ALPHA,PSI,LAMBDA,P,INVCOR,COMM,RESCOR,GSCALE, 0000811 + IWORK,WORK3,IEROFC(NF)) 0000812 LNLIKE=(WORK3(1)+LNDET+P+P*LN2PI)*DBLE(N) 0000813 DF=IWORK(1) 0000814 NITERS=IWORK(2) 0000815 FMIN=WORK3(1) 0000816 CHISQ=WORK3(3) 0000817 C IV=1 ALLOWS COMMUNALITIES TO BE REUSED 0000818 IV=1 0000819 DO 20,I=1,P 0000820 OLDPSI(I,NF)=PSI(I) 0000821 20 CONTINUE 0000822 RETURN 0000823 END 0000824 ************************************************************************0000825 C SUBROUTINE OFCOMM(R,NV,NF,IND,NT,IV,MAXIT,MAXTRY, 0000826 C + EPS,EPSE,ALPHA,V,A,IA,RI,Y,S,G,IS,WK,IER) 0000827 C IMSL ROUTINE FOR COMPUTING UNROTATED FACTOR LOADING MATRIX 0000828 C WARNING ERRORS: 0000829 C IER=39 INDICATES THAT THE NUMBER OF DEGREES OF FREEDOM OF THE 0000830 C CHI-SQUARED STATISTIC WAS NOT POSITIVE. THE TEST WAS 0000831 C SKIPPED. ANY RESULTS SHOULD BE SUSPECT. 0000832 C IER=40 INDICATES MDCH COULD NOT CALCULATE THE PROBABILITY ALPHA0000833 C ASSOCIATED WITH IS(1) DEGREES OF FREEDOM WITH VALUE 0000834 C WK(3). ONLY OUTPUT VARIABLE ALPHA IS AFFECTED. 0000835 C IER=41 INDICATES SIGNIFICANCE LEVEL NOT ATTAINED IN MODEL WITH 0000836 C NF FACTORS. 0000837 C IER=42 INDICATES THAT THE HEYWOOD CASE WAS ENCOUNTERED ON THE 0000838 C ITERATION THAT A TERMINAL ERROR WAS DETECTED. 0000839 C TERMINAL ERRORS: 0000840 C IER=129 INDICATES AT LEAST ONE OF NV, NF, IA, OR IND WAS 0000841 C SPECIFIED INCORRECTLY. 0000842 C IER=130 INDICATES R WAS NOT POSITIVE DEFINITE. 0000843 C IER=131 INDICATES MAXIT WAS EXCEEDED. 0000844 C IER=132 INDICATES MAXTRY WAS EXCEEDED. 0000845 C IER=133 INDICATES EIGRS DID NOT CONVERGE. CHECK DATA FOR 0000846 C REDUNDANCY OR INCREASE THE PRECISION TO DOUBLE. 0000847 C IER=134 INDICATES A SECOND DERIVATIVE MATRIX WAS NOT POSITIVE 0000848 C DEFINITE. CHECK THE DATA FOR REDUNDANCY AND CHECK THAT0000849 C THE NUMBER OF DEGREES OF FREEDOM IS(1) IS POSITIVE. 0000850 C END 0000851 ************************************************************************0000852 SUBROUTINE DOCRIT(N,P,NF,LNLIKE,BESTNF,DOAIC,NPARMS,PENALT,AIC, 0000853 + CAIC,CRITER) 0000854 INTEGER P,NF,BESTNF,N 0000855 DOUBLE PRECISION LNLIKE,NPARMS,PENALT,AIC,CAIC,CRITER(0:P,5) 0000856 LOGICAL DOAIC 0000857 C LOCAL VARIABLES 0000858 DOUBLE PRECISION TWO,MINAIC,MINCA,PENAIC,PENCA 0000859 PARAMETER(TWO=2.0D0) 0000860 INTRINSIC DBLE,LOG 0000861 SAVE MINAIC,MINCA 0000862 DATA MINAIC,MINCA/1.0D38,1.0D38/ 0000863 IF(NF.EQ.0) THEN 0000864 NPARMS=0.5D0*P*(P+1) 0000865 CRITER(0,1)=NPARMS 0000866 ELSE 0000867 NPARMS=DBLE(NF*P+P-0.5D0*NF*(NF-1.0D0)) 0000868 CRITER(NF,1)=NPARMS 0000869 ENDIF 0000870 PENAIC=TWO*NPARMS 0000871 AIC=LNLIKE+PENAIC 0000872 CRITER(NF,2)=AIC 0000873 IF (DOAIC) THEN 0000874 IF (AIC.LT.MINAIC) THEN 0000875 MINAIC=AIC 0000876 BESTNF=NF 0000877 ENDIF 0000878 PENALT=PENAIC 0000879 ENDIF 0000880 PENCA=LOG(DBLE(N))*NPARMS 0000881 CAIC=LNLIKE+PENCA 0000882 CRITER(NF,3)=CAIC 0000883 IF (.NOT.DOAIC) THEN 0000884 IF (CAIC.LT.MINCA) THEN 0000885 MINCA=CAIC 0000886 BESTNF=NF 0000887 ENDIF 0000888 PENALT=PENCA 0000889 ENDIF 0000890 RETURN 0000891 END 0000892 ************************************************************************0000893 SUBROUTINE ROTATE(P,NF,LAMBDA,WORK1,WORK2,WORK3) 0000894 INTEGER P,NF 0000895 DOUBLE PRECISION LAMBDA(P,P), 0000896 + WORK1(P,P*(P+1)/2), 0000897 + WORK3(P),WORK2(P,P) 0000898 C LOCAL VARIABLES 0000899 INTEGER IER,I,J,NORM,II,MAXIT,PLUS,NEG 0000900 DOUBLE PRECISION W,EPS,DELTA,ZERO 0000901 EXTERNAL OFROTA 0000902 C CONVERGENCE CRITERIA FOR THE IMSL ROUTINE OFROTA: 0000903 C MAXIT: MAXIMUM NUMBER OF ITERATIONS 0000904 C IMSL SUGGESTS 30 0000905 C EPS: CONVERGENCE CONSTANT FOR ROTATION 0000906 C IMSL SUGGESTS 0.0001 0000907 C DELTA: CONVERGENCE CONSTANT FOR ROTATION 0000908 C IMSL SUGGESTS 0.001 0000909 PARAMETER(MAXIT=60, 0000910 + EPS=0.0001D0, 0000911 + DELTA=0.001D0) 0000912 ZERO=0.0D0 0000913 C NORM=1 INDICATES THAT ROW NORMALIZATION IS PERFORMED IN OFROTA 0000914 NORM=1 0000915 C II=0 INDICATES AN IMAGE ANALYSIS IS NOT PERFORMED IN OFROTA 0000916 II=0 0000917 C W=1.0 INDICATES THAT THE VARIMAX METHOD IS USED IN OFROTA 0000918 W=1.0D0 0000919 C COMPUTING ORTHOGONAL ROTATION OF LAMBDA 0000920 CALL OFROTA(LAMBDA,P,P,NF,NORM,II,MAXIT,W,EPS,DELTA,LAMBDA,P, 0000921 + WORK2,P,WORK3,WORK1,IER) 0000922 C MAKE CORRELATIONS MOSTLY POSITIVE 0000923 DO 30, J=1,NF 0000924 PLUS=0 0000925 NEG=0 0000926 DO 10, I=1,P 0000927 IF (LAMBDA(I,J).LT.ZERO) THEN 0000928 NEG=NEG+1 0000929 ELSE 0000930 PLUS=PLUS+1 0000931 ENDIF 0000932 10 CONTINUE 0000933 IF (NEG.GT.PLUS) THEN 0000934 DO 20, I=1,P 0000935 LAMBDA(I,J)=-LAMBDA(I,J) 0000936 20 CONTINUE 0000937 ENDIF 0000938 30 CONTINUE 0000939 RETURN 0000940 C END 0000941 ************************************************************************0000942 C SUBROUTINE OFROTA(A,IA,NV,NF,NORM,II,MAXIT,W, 0000943 C + EPS,DELTA,B,IB,T,IT,F,WK,IER) 0000944 C IMSL ROUTINE FOR ORTHOGONAL ROTATION OF A FACTOR LOADING MATRIX 0000945 C TERMINAL ERROR: 0000946 C IER=129 INDICATES THAT AT LEAST ONE OF NV, NF, IA, IT, OR IB 0000947 C WAS SPECIFIED INCORRECTLY. 0000948 ************************************************************************0000949 END 0000950 ************************************************************************ APPENDIX E: THE CALL STATEMENTS FROM FACAIC 1: PROGRAM MAIN 62: CALL START(N,P,NROW,NCOL,IOUT,VERSON,DOAIC,RAWDAT,AUTO,SETNF) 64: CALL DESRD(N,P,X,MEAN,CORR,STD,IOUT,NROW,WORK2) 65: CALL DESOUT(MEAN,CORR,STD,P,IOUT) 67: CALL GETSYM(CORR,P,IOUT) 68: CALL SYMOUT(CORR,P,IOUT) 72: CALL FACAIC(N,P,DOAIC,CORR,AUTO,SETNF,BESTNF,ALPHA, 73: + NPARMS,PENALT,AIC,CAIC,LNDET,PSI,CRITER,BIGFAC, 74: + LAMBDA,INVCOR,COMM,RESCOR,IEROFC,LNLIKE,DF,NITERS, 75: + FMIN,CHISQ,OLDPSI,GSCALE,IWORK,WORK1,WORK2,WORK3) 76: CALL PRTOUT(N,P,COMM,LNDET,PSI,LAMBDA,BESTNF,RESCOR,IOUT, 77: + ALPHA,NPARMS,PENALT,AIC,CAIC,CRITER,BIGFAC,SETNF,AUTO, 78: + BESTNF,DOAIC,IEROFC,LNLIKE,DF,NITERS,FMIN,CHISQ, 79: + CORR,WORK1,WORK2) 86: SUBROUTINE START(N,P,NROW,NCOL,IOUT,VERSON,DOAIC,RAWDAT,AUTO, 87: + SETNF) 101: CALL GTOPFN(UN,FN,NEW,OLD,TITLE) 105: CALL UGETIO(IOPT,NIN,IOUT) 107: CALL UGETIO(IOPT,NIN,IOUT) 132: CALL QUERY(TITLE,DOAIC) 144: CALL GTOPFN(UN,FN,NEW,OLD,TITLE) 151: CALL QUERY(TITLE,RAWDAT) 153: CALL QUERY(TITLE,AUTO) 169: SUBROUTINE QUERY(TITLE,RESPON) 188: SUBROUTINE GTOPFN(UN,FN,NEW,OLD,TITLE) 234: SUBROUTINE GETSYM(X,P,IOUT) 246: CALL QUERY(TITLE,FREE) 268: SUBROUTINE SYMOUT(X,P,IOUT) 280: CALL USWSM(ITITLE,NC,X,P,IOPT) 285: SUBROUTINE DESRD(N,P,X,MEAN,CORR,STD,IOUT,NROW,WORK2) 298: CALL QUERY(TITLE,FREE) 313: CALL GETDAT(FREE,X,NROW,P,IN,INUNIT,IOUT) 319: CALL BECOVW(X,NROW,WT,NBR,XMISS,MEAN,CORR,INCD,WORK2,IER) 347: SUBROUTINE GETDAT(FREE,X,NROW,P,IN,INUNIT,IOUT) 388: SUBROUTINE DESOUT(MEAN,CORR,STD,P,IOUT) 401: CALL USWFV(ITITLE,NC,MEAN,P,INC,IOPT) 407: CALL USWFV(ITITLE,NC,STD,P,INC,IOPT) 412: CALL USWSM(ITITLE,NC,CORR,P,IOPT) 425: SUBROUTINE PRTOUT(N,P,COMM,LNDET,PSI,LAMBDA,NF,RESCOR,IOUT,ALPHA, 426: + NPARMS,PENALT,AIC,CAIC,CRITER,BIGFAC,SETNF,AUTO, 427: + BESTNF,DOAIC,IEROFC,LNLIKE,DF,NITERS,FMIN,CHISQ, 428: + CORR,WORK1,WORK2) 466: CALL USWFV(ITITLE,NC,COMM,P,INC,IOPT) 473: CALL USWFV(ITITLE,NC,PSI,P,INC,IOPT) 479: CALL USWFM(ITITLE,NC,LAMBDA,P,P,NF,IOPT) 482: CALL VMULFP(LAMBDA,LAMBDA,P,NF,P,P,P,WORK1,P,IER) 487: CALL VCVTFS(WORK1,P,P,WORK2) 493: CALL USWSM(ITITLE,NC,WORK2,P,IOPT) 499: CALL USWSM(ITITLE,NC,RESCOR,P,IOPT) 509: CALL USWSM(ITITLE,NC,WORK2,P,IOPT) 514: CALL LINV3F(WORK1,B,IJOB,P,P,D1,D2,WORK2,IER) 586: SUBROUTINE FACAIC(N,P,DOAIC,CORR,AUTO,SETNF,BESTNF,ALPHA, 587: + NPARMS,PENALT,AIC,CAIC,LNDET,PSI,CRITER,BIGFAC, 588: + LAMBDA,INVCOR,COMM,RESCOR,IEROFC,LNLIKE,DF, 589: + NITERS,FMIN,CHISQ,OLDPSI,GSCALE, 590: + IWORK,WORK1,WORK2,WORK3) 701: CALL UERSET(LEVEL,LEVOLD) 708: CALL LUDECP(CORR,WORK1,P,D1,D2,IER) 714: CALL DOCRIT(N,P,IZERO,ZERLIK,BESTNF,DOAIC,NPARMS,PENALT,AIC,CAIC, 715: + CRITER) 719: CALL FACTOR(N,P,NF,CORR,LNDET,PSI,LAMBDA,INVCOR,COMM, 720: + RESCOR,LNLIKE,DF,NITERS,FMIN,CHISQ,FIRST,ALPHA, 721: + IEROFC,GSCALE,OLDPSI,IWORK,WORK3) 726: CALL DOCRIT(N,P,NF,LNLIKE,BESTNF,DOAIC,NPARMS,PENALT,AIC, 727: + CAIC,CRITER) 733: CALL FACTOR(N,P,NF,CORR,LNDET,PSI,LAMBDA,INVCOR,COMM, 734: + RESCOR,LNLIKE,DF,NITERS,FMIN,CHISQ,FIRST,ALPHA, 735: + IEROFC,GSCALE,OLDPSI,IWORK,WORK3) 740: CALL DOCRIT(N,P,NF,LNLIKE,BESTNF,DOAIC,NPARMS,PENALT,AIC, 741: + CAIC,CRITER) 748: CALL FACTOR(N,P,BESTNF,CORR,LNDET,PSI,LAMBDA,INVCOR,COMM, 749: + RESCOR,LNLIKE,DF,NITERS,FMIN,CHISQ,FIRST,ALPHA, 750: + IEROFC,GSCALE,OLDPSI,IWORK,WORK3) 752: CALL DOCRIT(N,P,BESTNF,LNLIKE,BESTNF,DOAIC,NPARMS,PENALT,AIC, 753: + CAIC,CRITER) 755: CALL ROTATE(P,BESTNF,LAMBDA,WORK1,WORK2,WORK3) 771: SUBROUTINE FACTOR(N,P,NF,CORR,LNDET,PSI,LAMBDA,INVCOR,COMM, 772: + RESCOR,LNLIKE,DF,NITERS,FMIN,CHISQ,FIRST,ALPHA, 773: + IEROFC,GSCALE,OLDPSI,IWORK,WORK3) 809: CALL OFCOMM(CORR,P,NF,IND,N,IV,MAXIT,MAXTRY,EPS,EPSE, 810: + ALPHA,PSI,LAMBDA,P,INVCOR,COMM,RESCOR,GSCALE, 811: + IWORK,WORK3,IEROFC(NF)) 852: SUBROUTINE DOCRIT(N,P,NF,LNLIKE,BESTNF,DOAIC,NPARMS,PENALT,AIC, 853: + CAIC,CRITER) 893: SUBROUTINE ROTATE(P,NF,LAMBDA,WORK1,WORK2,WORK3) 920: CALL OFROTA(LAMBDA,P,P,NF,NORM,II,MAXIT,W,EPS,DELTA,LAMBDA,P, 921: + WORK2,P,WORK3,WORK1,IER) The following subroutines are from IMSL: UGETIO USWSM BECOVW USWFV USWFM VMULFP VCVTFS LINV3F UERSET LUDECP OFCOMM OFROTA ************************************************************************