C********** MSFORTRAN77 ADAPTATION OF CLOGG LATENT CLASS PROGRAM C LARGE MODEL C ADAPTED TO RUN ON IBM PC's AND COMPATABLES BY RANDALL KNACK C MARCH, 1986 C PROGRAM FOR CALCULATION OF MAXIMUM LIKELIHOOD LATENT CLASS MODEL C WILL HANDLE L. C. MODEL WITH NX LATENT CLASSES, NVAR MANIFEST C VARIABLES, WITH 'IV' CLASSES PER VARIABLE. ONLY DIMENSION STATE- C MENT MUST BE CHANGED TO HANDLE GENERAL CASE. C PRESENT PROGRAM IS CONSIDERABLE REVISION OF PGM USED BY GOODMAN C IN 1974 ARTICLES C VERSION 5-- AUGUST 1982 C Clogg's manual (July, 1977) notes the following: C NVAR = # of manifest variables C TSIZE = product of classes of variables (eg. IxJxKxL) C IMAX = number of classes of the variable with the most classes C NX = # of latent classes C DIMENSION FMT(20),IV(NVAR),TITLE(20),VNAME(NVAR), C XLABEL(NVAR,IMAX),TM(NVAR,IMAX),N(NVAR),R(TSIZE) C DIMENSION X(NX),FS(TSIZE,FF(TSIZE),P(TSIZE),Q(TSIZE,NX), C Y(NVAR,NX,IMAX),S(TSIZE,NX),PQ(TSIZE,NX),TT(IMAX),IX(NX), C IY(NVAR,NX,IMAX),XT(NX),YT(NVAR,NX,IMAX) INTEGER TSIZE REAL*8 MAXDEV,DEV,FN,T,SUM,SUMY,DSQRT CHARACTER FMT*20 DIMENSION IV(6),TITLE(20),VNAME(6),XLABEL(6,6),TM(6,6), + N(6),TT(6),IX(20),IY(6,20,6),MARK(60) REAL*8 X(20),FS(243),FF(243),P(243),Q(243,20),Y(6,20,6),G(10), + S(243,20),PQ(243,20), XT(20),YT(6,20,6),R(243),FT(243) REAL*8 A(243,60) DATA IX/20*0/,IY/720*0/,TM/36*0.0/,A/14580*0.D0/ *add MS-DOS filehandling routine CALL FILES WRITE(6,10000) 10000 FORMAT(//35X,53H ************************************************* +***) WRITE(6,10004) 10004 FORMAT(35X,53H * LATENT STRUCTURE ANALYSIS -- MLLSA (CLOGG, 1977) +*) WRITE(6,10008) 10008 FORMAT(35X,53H * ADAPTED TO THE PC (R.KNACK, 1986) +*) WRITE(6,10010) 10010 FORMAT( 35X,53H ************************************************* +***) 1 READ(8,115,END=900)TITLE 115 FORMAT(20A4) DO 1235 IJ=1,60 MARK(IJ)=0 1235 CONTINUE WRITE (6,116)TITLE 116 FORMAT(//,13H RUN NAME: ,1H ,20A4) READ (8,100,END=900)NVAR,NX,FN,MAXIT,MAXDEV,IND,ILDIS,IADD,IVLAB, *ILAB,IXRES,IVARES,IT,IRESID,JOB,IDENT,NG 100 FORMAT (2I2,F6.0,I6,F8.7,12I2) IF(NVAR.EQ.0)STOP IF(JOB.NE.0)GO TO 50 READ(8,1000)(IV(J),J=1,NVAR) 1000 FORMAT(40I2) C LABELS FOR VARIABLES AND CATEGORIES PER VARIABLE IF(IVLAB.EQ.0)GO TO 1111 READ(8,1110)(VNAME(I),I=1,NVAR) IF(ILAB.EQ.0)GO TO 1112 DO 1200 I=1,NVAR JJ=IV(I) 1200 READ(8,1110)(XLABEL(I,J),J=1,JJ) 1110 FORMAT (10A8) GO TO 1300 1111 DO 1311 I=1,NVAR 1311 VNAME(I)=I 1112 DO 1310 I=1,NVAR JJ=IV(I) DO 1310 J=1,JJ 1310 XLABEL(I,J)=J 1300 TSIZE=IV(1) DO 1001 J=2,NVAR 1001 TSIZE=TSIZE*IV(J) READ(8,'(A)')FMT READ(8,FMT)(FS(I),I=1,TSIZE) READ(8,102)(X(I),I=1,NX) 22 T=0.D0 DO 2 I=1,TSIZE 2 T=T+FS(I) WRITE(6,'(/A,T55,F8.0)')'TOTAL N OF INPUT DATA (card 7.) =',T WRITE(6,'( A,T55,F8.0)')'EXPECTED N OF INPUT DATA (card 2, cols. 5 +-10) =',FN * WRITE (6,103)T,FN * 103 FORMAT(/60H TOTAL OF INPUT (card 7.), AND EXPECTED (card 2, cols. * + 5-10),2F8.0) 102 FORMAT(10F8.7) FN=T WRITE (6,104)(FS(I),I=1,TSIZE) 104 FORMAT(16F8.0) WRITE(6,'(/A)')'INPUT X''S are start values for the latent class p +robabilities (card 8.)' WRITE (6,105)(X(I),I=1,NX) 105 FORMAT (14H INPUT X'S ==>,10F10.7/) WRITE(6,'(/A)')'INPUT Y''S are start values for the conditional pr +obailities (card 9.)' WRITE(6,106) 106 FORMAT (10H INPUT Y'S) DO 3 J=1,NVAR KK=IV(J) READ(8,1071)((Y(J,I,K),K=1,KK),I=1,NX) 1071 FORMAT(10F8.7) WRITE(6,107)((Y(J,I,K),K=1,KK),I=1,NX) 3 CONTINUE 107 FORMAT (8F10.7) C CALCULATE IMAX = MAX CLASS IV(I) IMAX=IV(1) DO 1002 J=2,NVAR IF(IV(J).LE.IMAX)GO TO 1002 IMAX=IV(J) 1002 CONTINUE C COMPUTE CHI-SQUARE FOR INDEPENDENCE MODEL--OPTIONAL IF(IND.EQ.0)GO TO 50 K=1 DO 52 INDEX=1,NVAR K=K*IV(INDEX) DO 53 I=1,TSIZE J=I-1 J=((J-(J/K)*K)*IV(INDEX))/K + 1 53 TM(INDEX,J)=TM(INDEX,J)+FS(I) 52 CONTINUE DO 54 I=1,TSIZE 54 FF(I)=1.D0 K=1 DO 55 INDEX=1,NVAR K=K*IV(INDEX) DO 56 I=1,TSIZE J=I-1 J=((J-(J/K)*K)*IV(INDEX))/K + 1 IF(INDEX.GT.1)GO TO 566 FF(I)=FF(I)*TM(INDEX,J) GO TO 56 566 FF(I)=FF(I)*TM(INDEX,J)/FN 56 CONTINUE 55 CONTINUE CHISQP=0.D0 CHISQL=0.D0 DO 57 I=1,TSIZE IF(FF(I).LE.0.D0)GO TO 57 CHISQP=CHISQP+((FS(I)-FF(I))**2)/FF(I) IF(FS(I).GT.0.D0)CHISQL=CHISQL+2.D0*FS(I)*DLOG(FS(I)/FF(I)) 57 CONTINUE WRITE(6,580) 580 FORMAT(//24H MARGINALS FOR IND MODEL) JJ=IV(NVAR) DO 583 J=1,JJ 583 G(J)=TM(NVAR,J)/FN DO 581 I=1,NVAR JJ=IV(I) WRITE(6,582)I,(TM(I,J),J=1,JJ) DO 581 J=1,JJ TM(I,J)=0.0 581 CONTINUE 582 FORMAT(I2,5X,10F8.0) WRITE(6,58)CHISQL,CHISQP 58 FORMAT(/32H CHI-SQUARE FOR IND MODEL LR= ,F10.3,5X,8H PRSN= ,F1 +0.3) 50 CONTINUE IF(MAXIT.EQ.0)MAXIT=500 IF(MAXDEV.EQ.0.D0)MAXDEV=5.0D-05 ITER=1 C RESTRICTIONS: C RESTRICTIONS ON X ARE DEFINED BY IXRES NE 0 AND RESTRICTIONS ARE C ENTERED IN IX(NX) C RESTRICTIONS ON YT ARE DEFINED BY IVARRES NE 0 AND RESTRICTIONS C ARE ENTERED IN IY(NVAR,NX,IV(NVAR)) C THE NUMBER '0' (OR BLANK) DEFINES A FREE PARAMETER. C THE NUMBER '1' INDICATES THAT INITIAL Y OR INITIAL X IS RESTORED C AFTER EACH ITERATION. C INTEGERS GREATER THAN '1' FILL OUT THE REMAINING ELEMENTS OF IY C MATRIX AND IX VECTOR. IF ANY SET OF INTEGERS (NOT '0' OR '1') C ARE EQUAL TO EACH OTHER, THEN THESE PROBABILITIES ARE CON- C STRAINED TO BE EQUAL BEGINNING (OR ENDING) EACH ITERATION C AFTER THE FIRST. C BOTH Y AND X ARE RESCALED AFTER EACH ITERATION TO ENSURE THAT C SUM YT OVER CLASSES OF MANIFEST VAR=1, AND SUM XT=1. RESCAL- C ING ONLY AFFECTS FREE PARAMETERS CORRESPONDING TO IY EQ '0' C AND IX EQ '0'. FOR MOST APPLICATIONS THERE MUST BE ONE C FREE PARAMETER WITHIN EACH FIXED SET OF NVAR =I, NX = J. C EXCEPTIONS TO ABOVE RESCALING MAY OCCUR WITH TWO EQUALITY C RESTRICTIONS WITHIN THE GIVEN SET OF I,J. THESE ARE C DENOTED BY ENTERING RESTRICTIONS NUMBERS IN IY OF RANGE 30-39 C THE PROGRAM WILL RESCALE THESE VALUES SO THAT EQUALITY C CONSTRAINTS ARE MAINTAINED. IF MORE GENERAL KINDS OF CONSTRAINTS C ARE NECESSARY THEN THE PROGRAM SECTION ON EQUALITY CONSTRAINTS C MUST BE MODIFIED. C C ENTER X RESTRICTION S IF IXRES NE 0 WRITE(6,*) CHAR(10),'X RESTRICTIONS are restrictions on the latent + class probabilities (card 10.)' IF(IXRES.EQ.O) THEN WRITE(6,*)'X RESTRICTIONS ==> NONE' ELSE WRITE(6,125) 125 FORMAT(15H X RESTRICTIONS) READ(8,71)(IX(I),I=1,NX) WRITE(6,701)(IX(I),I=1,NX) 71 FORMAT(20I2) 701 FORMAT(20I4) ENDIF C C ENTER Y RESTRICTIONS IF IVARRES NE 0. NOTE THAT WE MUST ENTER C IV(I)*NX RESTRICTIONS FOR EACH VARIABLE I) WRITE(6,*) CHAR(10),'Y RESTRICTIONS are restrictions on the condit +ional probabilities (card 11.)' 722 WRITE(6,1125) 1125 FORMAT(15H Y RESTRICTIONS) IF(IVARES.EQ.0)GO TO 4 DO 73 I=1,NVAR KK=IV(I) 73 READ(8,74)((IY(I,J,K),K=1,KK),J=1,NX) 74 FORMAT(20I4) WRITE (6,1127) 1127 FORMAT(17H LATENT CLASS ==>,5X,2H 1,2X,2H 2,2X,2H 3,2X,2H 4) WRITE(6,*)'ITEM # RESPONSE' DO 705 I=1,NVAR KK=IV(I) DO 705 K=1,KK IF((IVLAB.NE.0).AND.(ILAB.NE.0))WRITE(6,1126)VNAME(I),XLABEL(I,K), *(IY(I,J,K),J=1,NX) IF(IVLAB.EQ.0)WRITE(6,1128)VNAME(I),XLABEL(I,K),(IY(I,J,K),J=1,NX) IF((IVLAB.NE.0).AND.(ILAB.EQ.0))WRITE(6,1129)VNAME(I),XLABEL(I,K), *(IY(I,J,K),J=1,NX) 1126 FORMAT(A8,2X,A8,2X,10I4,/,20X,10I4,/,20X,10I4) 1128 FORMAT(F3.0,7X,F3.0,7X,10I4,/,20X,10I4,/,20X,10I4) 1129 FORMAT(A8,2X,F3.0,7X,10I4,/,20X,10I4,/,20X,10I4) 705 CONTINUE C C BEGIN CALCULATION OF PARAMETERS 4 DO 5 I=1,NVAR 5 N(I)=1 DO 10 INDEX=1,TSIZE DO 6 J=1,NX Q(INDEX,J)=1. DO 6 I=1,NVAR Q(INDEX,J)=Q(INDEX,J)*Y(I,J,N(I)) 6 CONTINUE I=1 7 N(I)=N(I)+1 IF(N(I).LE.IV(I))GO TO 10 N(I)=1 I=I+1 IF(I.LE.NVAR)GO TO 7 10 CONTINUE DO 13 INDEX=1,TSIZE P(INDEX)=0.D0 DO 11 I=1,NX IF(Q(INDEX,I).GT.1.D-12)PQ(INDEX,I)=Q(INDEX,I)*X(I) IF(Q(INDEX,I).LE.1.D-12)PQ(INDEX,I)=0.D0 IF((PQ(INDEX,I).LE.0.D0).OR.(Q(INDEX,I).LE.1.D-07))PQ(INDEX,I)= 1 0.D0 11 P(INDEX)=P(INDEX)+PQ(INDEX,I) DO 12 I=1,NX S(INDEX,I)=0.D0 IF(P(INDEX).LE..0000001)GO TO 12 S(INDEX,I)=PQ(INDEX,I)/P(INDEX) 12 CONTINUE 13 FF(INDEX)=FN*P(INDEX) IF(ITER.GT.1.OR.IT.EQ.0)GO TO 15 WRITE(6,108) 108 FORMAT(///32H FIRST P'S=PROB OF MANIFEST CELL) WRITE(6,107)(P(I),I=1,TSIZE) WRITE(6,7107) 7107 FORMAT(/67H FIRST PQ'S=JOINT PROBABILITY OF MANIFEST CELL I AND LA +TENT CLASS J) DO 75 I=1,NX WRITE(6,7109)I 7109 FORMAT(14H LATENT CLASS=,I4) 75 WRITE(6,7108)(PQ(J,I),J=1,TSIZE) 7108 FORMAT(8X,8F10.7) WRITE(6,171) 171 FORMAT(/10H FIRST Q'S) DO 76 I=1,NX 76 WRITE(6,107)(Q(J,I),J=1,TSIZE) WRITE(6,1077) 1077 FORMAT(/10H FIRST S'S) DO 14 I=1,NX 14 WRITE(6,107)(S(J,I),J=1,TSIZE) WRITE(6,109) 109 FORMAT(/16H FIRST ESTIMATES) WRITE(6,110)(FF(I),I=1,TSIZE) 110 FORMAT(8F10.2) ISTART=1 GO TO 60 19 WRITE(6,130)CHISQL,CHISQP 130 FORMAT(/11H INITIAL LR,F10.2,5X,6H PRSN ,F10.2) 15 CONTINUE DEV=0.D0 DO 16 INDEX=1,NX T=0.D0 DO 17 I=1,TSIZE 17 T=T+FS(I)*S(I,INDEX) XT(INDEX)=T/FN C C INSERT TEST OF LATENT CLASS PROPORTION APPROACHING ZERO C PROGRAM WILL STOP PRINTING A WARNING AND THE LATENT CLASS C WHERE THE PROBLEM AROSE. IF(XT(INDEX).GE..0000001) GO TO 16 MAXIT=0 WRITE(6,7110)ITER,INDEX 7110 FORMAT (/54H WARNING--LATENT CLASS PROPORTION WAS LESS THAN 1.0E-7, * 5X,11HITERATION= ,I4,5X,21HLATENT CLASS NUMBER= ,I4) C 16 CONTINUE K=1 DO 25 INDEX=1,NVAR K=K*IV(INDEX) DO 23 M=1,NX DO 2202 L=1,IMAX 2202 TT(L)=0.D0 DO 21 I=1,TSIZE J=I-1 J=((J-(J/K)*K)*IV(INDEX))/K + 1 TT(J)=TT(J)+FS(I)*S(I,M) 21 CONTINUE C PSSIBLE TROUBLE SPOT LL=IV(INDEX) DO 222 L=1,LL 222 YT(INDEX,M,L)=TT(L)/(FN*XT(M)) 23 CONTINUE 25 CONTINUE C C TEST FOR Y RESTRICTIONS C IF(IVARES.EQ.0)GO TO 30 DO 31 I=1,NVAR DO 31 J=1,NX KK=IV(I) DO 31 K=1,KK IF(IY(I,J,K).EQ.1)YT(I,J,K)=Y(I,J,K) IF(IY(I,J,K).LE.1)GO TO 31 SUM=0.D0 T=0.D0 IP=I II=IY(I,J,K) DO 36 L=IP,NVAR KKK=IV(L) DO 360 M=1,NX DO 361 IN=1,KKK IF(II.NE.IY(L,M,IN))GO TO 361 SUM=SUM+YT(L,M,IN)*XT(M) T=T+XT(M) IY(L,M,IN)=-IY(L,M,IN) 361 CONTINUE 360 CONTINUE 36 CONTINUE II=IY(I,J,K) DO 37 LL=IP,NVAR KKK=IV(LL) DO 370 MM=1,NX DO 371 IIN=1,KKK IF(II.NE.IY(LL,MM,IIN))GO TO 371 YT(LL,MM,IIN)=SUM/T 371 CONTINUE 370 CONTINUE 37 CONTINUE 31 CONTINUE C C RESCALE Y SO THAT SUM YT(I,J,K) OVER K = 1.0 C RESTORE MATRIX OF RESTRICTIONS C NOTE THAT ONE 'FREE' PARAMETER IS NECESSARY FOR COND'L PROBS C WHERE LATENT CLASS IS FIXED. NECESSARY TO RESCALING 30 DO 389 I=1,NVAR DO 389 J=1,NX SUM=0.D0 SUMY=0.D0 KK=IV(I) DO 38 K=1,KK IY(I,J,K)=IABS(IY(I,J,K)) SUM=SUM+YT(I,J,K) IF (IY(I,J,K).GT.0)SUMY=SUMY+YT(I,J,K) 38 CONTINUE C IF NO FREE PARMS, THEN SUM-SUMY = 0.0 AND (IF RESTRICTIONS ARE NOT C REDUNDANT) THERE MUST BE EQUALITY RESTRICTION WITHIN I'TH VAR C AND J'TH LATENT CLASS. THIS PART WILL RESCALE OK WHEN TWO ARE C EQUAL, BUT NEEDS MODIFICATION IF MORE THAN TWO ARE EQUAL. C A SPECIAL SIGNAL OF TWO DIGIT NUMBER BEGINNING WITH 3 KEYS C THE POSITION OF THE WITHIN (I,J) EQUALITY RESTRICTIONS. IF ((SUM-SUMY))381,381,382 381 DO 383 K=1,KK IF((IY(I,J,K)/30).EQ.1)YT(I,J,K)=.5 - SUM/2.0 - YT(I,J,K) 383 CONTINUE GO TO 389 382 DO 388 K=1,KK IF (IY(I,J,K).EQ.0)YT(I,J,K)=YT(I,J,K)*(1.-SUMY)/(SUM-SUMY) 388 CONTINUE 389 CONTINUE C C TEST FOR X RESTRICTIONS IF(IXRES.EQ.0)GO TO 810 IF(NG.NE.0)GO TO 801 DO 744 L=1,NX IF(IX(L).EQ.1)XT(L)=X(L) IF((IX(L).LE.1).OR.(L.EQ.NX))GO TO 744 SUM=XT(L) K=L+1 XK=1. DO 766 I=K,NX IF(IX(L).NE.IX(I))GO TO 766 XK=XK+1. SUM=SUM+XT(I) 766 CONTINUE XT(L)=SUM/XK DO 77 I=K,NX IF(IX(L).NE.IX(I))GO TO 77 IX(I)=-IX(I) XT(I)=SUM/XK 77 CONTINUE 744 CONTINUE GO TO 810 801 CONTINUE DO 802 L=1,NX IF(IX(L).EQ.1)XT(L)=X(L) IF((IX(L).LE.1).OR.(L.EQ.NX))GO TO 802 NGG=(L-1)/(NX/NG) + 1 SUM=XT(L)/G(NGG) K=L+1 XK=1.D0 DO 804 I=K,NX IF(IX(L).NE.IX(I))GO TO 804 NGGG=(I-1)/(NX/NG) + 1 XK=XK + 1.D0 SUM= SUM + XT(I)/G(NGGG) 804 CONTINUE XT(L)=(SUM/XK)*G(NGG) DO 805 I=K,NX IF(IX(L).NE.IX(I))GO TO 805 NGG=(I-1)/(NX/NG) + 1 IX(I)=-IX(I) XT(I)=(SUM/XK)*G(NGG) 805 CONTINUE 802 CONTINUE C NOW RESCALE XT(I) SO THEY SUM TO 1.0 C RESTORE VECTOR OF X RESTRICTIONS JJ=1 DO 806 I=1,NG SUM=0.D0 SUMY=0.D0 JX=I*(NX/NG) DO 807 J=JJ,JX IX(J)=IABS(IX(J)) SUM=SUM + XT(J) IF(IX(J).NE.0)GO TO 807 SUMY=SUMY + XT(J) 807 CONTINUE DO 808 J=JJ,JX IF(IX(J).NE.0)GO TO 808 XT(J)=XT(J)*(G(I) - SUM + SUMY)/SUMY 808 CONTINUE 806 JJ=JX+1 GO TO 809 810 SUM=0.D0 SUMY=0.D0 DO 707 I=1,NX IX(I)=IABS(IX(I)) SUM=SUM+XT(I) IF(IX(I).NE.1)GO TO 707 SUMY=SUMY+XT(I) 707 CONTINUE DO 79 I=1,NX IF(IX(I).EQ.1)GO TO 79 XT(I)=XT(I)*(1.-SUMY)/(SUM-SUMY) 79 CONTINUE C NOW RESTORE VALUES TO BEGIN ITERATION, AND CHECK DEV 809 CONTINUE DO 92 K=1,NX IF(DABS(X(K)-XT(K)).GT.DEV)DEV=DABS(X(K)-XT(K)) X(K)=XT(K) DO 92 I=1,NVAR JJ=IV(I) DO 92 J=1,JJ IF(DABS(Y(I,K,J)-YT(I,K,J)).GT.DEV)DEV=DABS(Y(I,K,J)-YT(I,K,J)) 92 Y(I,K,J)=YT(I,K,J) IF((ITER/20)*20.EQ.ITER)THEN WRITE(6,112)DEV,ITER ENDIF IF(DEV.LE.MAXDEV)GO TO 41 IF(ITER.GT.MAXIT)GO TO 40 ITER=ITER+1 GO TO 4 40 WRITE(6,111) 111 FORMAT(/15H NO CONVERGENCE) 41 WRITE(6,112)DEV,ITER 112 FORMAT(/11H DEVIATION ,E16.8,5X,12H ITERATIONS ,I5) WRITE (6,113) 113 FORMAT(/4H FIT) WRITE(6,110)(FF(I),I=1,TSIZE) ISTART=2 60 CHISQL=0.D0 CHISQP=0.D0 DELTA=0.D0 IF(NG.NE.0)NGG=TSIZE/NG DO 61 I=1,TSIZE DELTA=DELTA+DABS(FF(I)-FS(I)) IF(IRESID.NE.0.AND.FF(I).GT.0.D0)R(I)=(FS(I)-FF(I))/DSQRT(FF(I)) FT(I)=DSQRT(FS(I))+DSQRT(FS(I)+1)-DSQRT(4.D0*FF(I)+1) IF(FF(I).LE.0.D0)GO TO 61 CHISQP=CHISQP+(FF(I)-FS(I))**2/FF(I) IF(FS(I).GT.0.D0)CHISQL=CHISQL+2.D0*FS(I)*DLOG(FS(I)/FF(I)) IF(NG.EQ.0)GO TO 61 IF((I/NGG)*NGG.EQ.I)THEN WRITE(6,121)CHISQL,CHISQP 121 FORMAT(//2H0 ,' CUMULATED CHISQUARE, LR=',F10.3,' PEARSON=',F10. +3) ENDIF 61 CONTINUE DELTA=DELTA/(2.D0*FN) IF(IRESID.EQ.0)GO TO 6616 WRITE(6,1130) 1130 FORMAT(//42H CELL OBSERVED EXPECTED (STDIZED RESID),16H (FREEMA +N-TUKEY)) DO 1132 I=1,TSIZE 1132 WRITE (6,1131)I,FS(I),FF(I),R(I),FT(I) 1131 FORMAT(I4,2X,F8.0,1X,F10.2,2H (,F6.2,2H ),6X,2H (,F6.2,2H )) 6616 IF(ISTART.EQ.1)GO TO 19 C WRITE(6,120)CHISQL,CHISQP,DELTA 120 FORMAT(/9H FINAL LR,F10.3,5X,6H PRSN ,F10.3,5X,23H INDEX OF DISSIM +ILARITY,F10.4) WRITE(6,300) 300 FORMAT(/33H FINAL LATENT CLASS PROBABILITIES) WRITE(6,107)(X(I),I=1,NX) WRITE(6,301) 301 FORMAT(/32H FINAL CONDITIONAL PROBABILITIES) WRITE(6,602) 602 FORMAT(17H LATENT CLASS ==>,9X,2H 1,6X,2H 2,6X,2H 3,6X,2H 4) DO 62 I=1,NVAR KK=IV(I) DO 62 K=1,KK IF((IVLAB.NE.0).AND.(ILAB.NE.0))WRITE(6,603)VNAME(I),XLABEL(I,K), +(Y(I,J,K),J=1,NX) IF(IVLAB.EQ.0)WRITE(6,604)VNAME(I),XLABEL(I,K),(Y(I,J,K),J=1,NX) IF((IVLAB.NE.0).AND.(ILAB.EQ.0))WRITE(6,605)VNAME(I),XLABEL(I,K), *(Y(I,J,K),J=1,NX) 603 FORMAT(A8,2X,A8,4X,10F8.4,/,22X,10F8.4,/,22X,10F8.4) 604 FORMAT(F3.0,7X,F3.0,10X,10F8.4,/,22X,10F8.4,/,22X,10F8.4) 605 FORMAT(A8,2X,F3.0,9X,10F8.4,/,22X,10F8.4,/,22X,10F8.4) 62 CONTINUE IF(IT.EQ.0)GO TO 909 WRITE(6,122) 122 FORMAT(/51H FINAL P'S=ESTIMATED PROBABILITY OF MANIFEST CELL I) WRITE(6,107)(P(I),I=1,TSIZE) WRITE(6,304) 304 FORMAT(/67H FINAL PQ'S=JOINT PROBABILITY OF MANIFEST CELL I AND LA +TENT CLASS J) DO 733 J=1,NX WRITE(6,7109)J 733 WRITE(6,7108)(PQ(I,J),I=1,TSIZE) WRITE(6,305) 305 FORMAT(/62H FINAL Q'S=PROBABILITY OF MANIFEST CELL I GIVEN LATENT +CLASS J) DO 63 J=1,NX WRITE(6,7109)J 63 WRITE(6,7108)(Q(I,J),I=1,TSIZE) WRITE(6,306) 306 FORMAT(/62H FINAL S'S=PROBABILITY OF LATENT CLASS J GIVEN MANIFEST + CELL I) DO 64 J=1,NX WRITE(6,7109)J 64 WRITE(6,7108)(S(I,J),I=1,TSIZE) C AUTOMATIC ASSIGNMENT OF RESPONDENTS TO LATENT CLASSES ON BASIS OF C MANIFEST RESPONSE--OPTIONAL, DEPENDING UPON ILDIS NE '0' 909 IF(ILDIS.EQ.0)GO TO 901 WRITE(6,902) 902 FORMAT(/70H ASSIGNMENT OF RESPONDENTS TO LATENT CLASS J GIVEN MANI +FEST RESPONSE I) WRITE(6,905) 905 FORMAT(55H CELL= OBSERVED= EXPECTED= ASSIGN TO CLASS= MODAL P=) XNUM=0.D0 DO 903 I=1,TSIZE MOD=1 CLMOD=S(I,1) DO 904 J=2,NX IF(S(I,J).LE.CLMOD)GO TO 904 MOD=J CLMOD=S(I,J) 904 CONTINUE XNUM=XNUM + FS(I)*CLMOD 903 WRITE(6,906)I,FS(I),FF(I),MOD,CLMOD PCT=(XNUM/FN)*100.D0 906 FORMAT(3X,I2,4X,F8.0,2X,F10.2,12X,I4,4X,F6.4) C COMPUTE 'LAMBDA' MEASURE OF ASSOCIATION BETWEEN X AND MANIF VARS CMODX=X(1) DO 908 I=2,NX IF(X(I).LE.CMODX)GO TO 908 CMODX=X(I) 908 CONTINUE XLMBDA=1.-(1.-PCT/100.)/(1.-CMODX) WRITE(6,907)PCT,XNUM,XLMBDA 907 FORMAT(/29H PERCENT CORRECTLY ALLOCATED=,F8.2,5X,28H NUMBER CORREC +TLY ALLOCATED=,F9.2,5X,8H LAMBDA=,F8.4) 901 IF(IDENT.EQ.0)GO TO 910 C C CALCULATION OF IDENTIFICATION OF MODEL PARAMETERS C 916 FORMAT(//,32H NUMBER OF ESTIMATED PARAMETERS=,I6,5X,36H DEGREES +OF FREEDOM IF IDENTIFIED=,I6) 918 FORMAT(/,'NUMBER OF ESTIMATED PARAMETERS GREATER THAN 60--CANNOT 1 USE SUBROUTINE GRAM-SCHMIDT TO CHECK IDENTIFICATION') IC=0 IROW=TSIZE-1 C BEGIN FILLING A MATRIX FOR DGRAM SUBROUTINE C FOR NOW, ONLY DEAL WITH CASE WHERE ALL PARAMETERES ARE EITHER C FREE OR FIXED. DEAL WITH EQUALITY RESTRICTIONS LATER C FIND FREE XT IM=0 DO 919 J=1,NX IF(IX(J).EQ.0)IM=J 919 CONTINUE DO 920 J=1,NX IF(XT(J).EQ.1.D0)GO TO 920 IF(J.EQ.IM)GO TO 920 IF(IX(J).EQ.1.OR.IX(J).LT.0)GO TO 920 IC=IC+1 DO 921 I=1,IROW 921 A(I,IC)=Q(I,J)-Q(I,IM) IF(IX(J).LE.1)GO TO 920 MARK(IC)=IX(J) IX(J)=-IX(J) 920 CONTINUE KK=1 DO 930 I=1,NVAR KK=KK*IV(I) DO 930 J=1,NX IP=IV(I) IM=0 DO 9301 K=1,IP IF(IY(I,J,K).EQ.0.AND.(YT(I,J,K).GT..5D-4.AND.YT(I,J,K).LT..99995) 1)IM=K 9301 CONTINUE DO 930 K=1,IP IF(IY(I,J,K).EQ.1.OR.IY(I,J,K).LT.0)GO TO 930 IF(YT(I,J,K).LT..5D-4.OR.YT(I,J,K).GT..99995)GO TO 930 IF(K.EQ.IM)GO TO 930 IC=IC+1 DO 931 M=1,IROW L=M-1 L=((L-(L/KK)*KK)*IV(I))/KK+1 IF(L.EQ.K)A(M,IC)=XT(J)*Q(M,J)/YT(I,J,L) IF(L.EQ.IM)A(M,IC)=-XT(J)*Q(M,J)/YT(I,J,L) 931 CONTINUE IF(IY(I,J,K).LE.1)GO TO 930 MARK(IC)=IY(I,J,K) IY(I,J,K)=-IY(I,J,K) 930 CONTINUE C FILL COLUMNS APPROPRIATELY, TAKING ACCOUNT OF EQUALITY CONSTRAINTS IC1=IC-1 DO 932 J=1,IC1 IF(MARK(J).LE.0)GO TO 932 J1=J+1 DO 933 JJ=J1,IC IF(MARK(JJ).NE.MARK(J))GO TO 933 MARK(JJ)=-1 DO 934 I=1,IROW 934 A(I,J)=A(I,J)+A(I,JJ) 933 CONTINUE 932 CONTINUE C NEG=0 DO 935 J=1,IC IF(MARK(J).EQ.-1)NEG=NEG+1 IF(MARK(J).EQ.-1)GO TO 935 MARK(J)=J-NEG 935 CONTINUE DO 936 J=1,IC IF(MARK(J).EQ.-1)GO TO 936 DO 937 I=1,IROW 937 A(I,MARK(J))=A(I,J) 936 CONTINUE IC=IC-NEG NDF=IROW-IC WRITE(6,916)IC,NDF C IF NUMBER OF ESTIMATED PARMS IS GT 60, CANNOT USE DGRAM IF(IC.LE.60)GO TO 917 WRITE(6,918) GO TO 910 917 IF(NDF.LT.0)GO TO 910 C C A MATRIX IS NOW FULL. WRITE IT OUT IF IT NE ZERO IF(IT.EQ.0)GO TO 943 DO 940 I=1,IROW 940 WRITE(6,942)(A(I,J),J=1,IC) 943 CALL DGRAM (A,243,IROW,IC,NR) NDF=IROW-NR WRITE(6,941)NR,NDF 941 FORMAT(//,13H COLUMN RANK=,I6,5X,20H DEGREES OF FREEDOM=,I6) 942 FORMAT(8F11.7) DO 944 I=1,NROW DO 944 J=1,IC 944 A(I,J)=0.D0 910 GO TO 1 900 STOP END C******** C******** C********** MSFORTRAN77 ADAPTATION OF CLOGG LATENT CLASS PROGRAM *$DEBUG $TITLE:'PC-Clogg BY R. Knack' $SUBTITLE:'SUBROUTINES' C ADAPTED TO RUN ON IBM PC's AND COMPATABLES BY RANDALL KNACK C MARCH, 1986 SUBROUTINE DGRAM(A,MR,NR,NC,IRANK) C C USING THE GRAM-SCHMIDT PROCESS, THIS SUBROUTINE COMPUTES C MUTUALLY ORTHOGONAL UNIT VECTORS WHICH SPAN THE SAME LINEAR C VECTOR SPACE AS THE INPUT MATRIX. IT IS BASED ON THE ALGORITHM C PRESENTED BY B. RUST IN CACM, MAY 1966, 9(5): 381-388. C C A IS THE MATRIX WHOSE VECTORS ARE TO BE ORTHOGONALIZED. C MR IS THE FIRST DIMENSION SUBSCRIPT OF A. C NR IS THE NUMBER OF ROWS IN A. C NC IS THE NUMBER OF COLUMNS IN A. C IRANK IS RETURNED AS THE RANK OF A. C DIMENSION A(MR,NC), U(60,60), AFLAG(60), ATEMP(60) REAL * 8 A,U,AFLAG,ATEMP,FAC,DOT,DOT1,DOT2,TOL,DMXSP,SQRT,DOTTOL DATA N /12/ SQRT(FAC) = DSQRT(FAC) DO 10 I = 1,NC DO 5 J = 1,NC 5 U(I,J) = 0. 10 U(I,I) = 1. C C NORMALIZE THE FIRST COL. OF A. C FAC = DMXSP(A(1,1),A(1,1),NR) *MSFORTRAN STOPS IF DIVISION BY ZERO NOT TRAPPED IF (FAC .EQ. 0.D00) GOTO 13 FAC = 1./SQRT(FAC) 13 DO 15 I = 1,NR 15 A(I,1) = A(I,1) * FAC DO 20 I = 1,NC 20 U(I,1) = U(I,1) * FAC AFLAG(1) = 1. C C OS/360 SINGLE PRECISION N=24; DOUBLE PRECISION N=56. C N IS THE DEPENDENT COL. TOL. FOR N BIT FLOATING POINT FRACTION. C TOL = (10. * .5**N)**2 C C PERFORM GRAMM-SCHMIDT ORTHOGONALIZATION PROCESS. C DO 100 J = 2,NC DOT1 = DMXSP(A(1,J),A(1,J),NR) JM1 = J - 1 DO 50 L = 1,2 DO 30 K = 1,JM1 30 ATEMP(K) = DMXSP(A(1,J),A(1,K),NR) DO 45 K = 1,JM1 DO 35 I = 1,NR 35 A(I,J) = A(I,J) - ATEMP(K) * A(I,K) * AFLAG(K) DO 40 I = 1,NC 40 U(I,J) = U(I,J) - ATEMP(K) * U(I,K) 45 CONTINUE 50 CONTINUE DOT2 = DMXSP(A(1,J),A(1,J),NR) *MSFORTRAN STOPS IF DIVISION BY ZERO NOT TRAPPED IF (DOT1 .EQ. 0.D00) THEN DOTTOL=0.D00 ELSE DOTTOL=DOT2/DOT1-TOL ENDIF IF(DOTTOL) 55,55,70 *OLD LINE * IF(DOT2/DOT1 - TOL) 55, 55, 70 55 DO 60 I = 1,JM1 ATEMP(I) = 0. DO 60 K = 1,I 60 ATEMP(I) = ATEMP(I) + U(K,I) * U(K,J) DO 65 I = 1,NR A(I,J) = 0. DO 65 K = 1,JM1 65 A(I,J) = A(I,J) - A(I,K) * ATEMP(K) * AFLAG(K) AFLAG(J) = 0. FAC = DMXSP(U(1,J),U(1,J),NC) FAC = 1./SQRT(FAC) GO TO 75 70 AFLAG(J) = 1. FAC = 1./SQRT(DOT2) 75 DO 80 I = 1,NR 80 A(I,J) = A(I,J) * FAC 100 CONTINUE IRANK=0 DO 110 J=1,NC IF(AFLAG(J).GT. .5D0) IRANK=IRANK+1 DO 110 I=1,NR 110 A(I,J)=A(I,J)*AFLAG(J) RETURN END C **************** C DMXSP C **************** DOUBLE PRECISION FUNCTION DMXSP(VECT1,VECT2,NELEM) C C DMXSP COMPUTES THE VECTOR PRODUCT OF THE VECTORS VECT1 AND VECT2 , C EACH OF WHICH HAS 'NELEM' ELEMENTS. C REAL*8 VECT1(NELEM),VECT2(NELEM) DMXSP=0.D00 DO 1000 I=1,NELEM 1000 DMXSP = DMXSP + VECT1(I)*VECT2(I) RETURN END C********************* C********************* ******************************************************************* *Randall E. Knack March, 1986 *Brief introduction and...file handling requests for MS-DOS. SUBROUTINE FILES IMPLICIT LOGICAL (A-Z) CHARACTER FNAME1*70, FNAME2*70, R*1, A(1)*80 LOGICAL*2 QEXIST A(1) = '========================================================== +=================' WRITE (*,*) A(1) WRITE (*,*) 'Unrestricted and Restricted Maximum Likelihood Latent + Structure Analysis' WRITE (*,*) 'An adaptation of MLLSA (Clifford Clogg) to the PC by +Randall Knack.' WRITE (*,*) 'Program setup follows Clogg ...Manual for Users (19 +77).' WRITE (*,*)'CURRENT PROGRAM LIMITS ARE AS FOLLOWS: (no more than)' WRITE (*,*)' 6 Manifest Variables' WRITE (*,*)' 20 Latent Classes' WRITE (*,*)' 6 Classes per Latent Variable' WRITE (*,*)' 243 Product of the number of classes for each la +tent variable' WRITE (*,*)' (ie. IxJxKxL=product given I...L = the numbe +r of classes for each' WRITE (*,*)' latent variable)' WRITE (*,*) A(1) 1000 WRITE (*,*) CHAR(10),'ENTER INPUT FILE NAME ==> ' READ (*,'(BN,A)',ERR=1000) FNAME1 INQUIRE (FILE=FNAME1,EXIST=QEXIST) IF (QEXIST) THEN OPEN (8,FILE=FNAME1,STATUS='OLD') ELSE GOTO 1000 ENDIF 1200 WRITE (*,*) CHAR(10),'ENTER OUTPUT FILE NAME ==> ' READ (*,'(BN,A)',ERR=1000) FNAME2 INQUIRE (FILE=FNAME2,EXIST=QEXIST) IF (QEXIST) THEN 1300 WRITE (*,*) 'FILE ALREADY EXISTS. OVERWRITE (Y/N) ?' READ (*,'(BN,A)',ERR=1300) R IF (R .EQ. 'Y' .OR. R .EQ. 'y') THEN OPEN (6,FILE=FNAME2,STATUS='OLD') ELSE GOTO 1200 ENDIF ELSE OPEN (6,FILE=FNAME2,STATUS='NEW') ENDIF RETURN END C****************************