C *** STAT PACK *** C SUBROUTINE FOR CHI SQUARES C CALLING SEQUENCE: CALL CHI(NV,NC,MV,MC,DATA,ISUB,LEVEL,NAMES) C WHERE NV - NUMBER OF VARIABLES ACTUALLY USED C NC - NUMBER OF OBSERVATION ACTUALLY USED C MV - MAX. NUMBER OF VARIABLES POSS. (SPEC. IN MAIN) C MC - MAX. NUMBER OF OBSERVATIONS POSS. (SPEC. IN MAIN) C DATA - STORAGE FOR DATA DIMENSIONED FOR MAXIMUM C ISUB IS A VECTOR DIMENSIONED AT LEAST FOR NC C LEVEL IS A VECTOR DIMENSIONED AT LEAST FOR NC C NAMES - IS A VECTOR CONTAINING VARIABLE NAMES C C CHI SQUARE OPTIONS ARE TO CALAPSE THE CONTINGENCY TABLE, AND C TO DETERMINE CERTAIN GROUPINGS OF CELLS PREVIOUS TO COLAPSEING C THE ROUTINE MAKES USE OF CHITAB, CHISOR, AND FISHER. C SUBROUTINE CHI(NV,NC,MV,MC,DATA,ISUB,LEVEL,NAMES) DIMENSION DATA(MC,MV),ISUB(1),LEVEL(1),IV(2),ANS(7),BRK(2,3,15) DIMENSION LVL(2),MCC(3,3),IVB(2) DIMENSION NAMES(1),STRG(50),B(5),LTBM(20) COMMON/DEV/ICC,IDATA,IOUT,IDLG,IDSK COMMON /PRNT/ LINPP,ICOPS,RUNPRG COMMON/EXTRA/HEDR(70),NSZ ILMAP="007777700000 IRMAP="000000077777 IJUST="000000100000 1 BREAK=0 CONTG=0 COLPS=0 FISHR=0 LVL(2)=0 IF(ICC.NE.2) WRITE(IDLG,3) 3 FORMAT(' LIST OPTIONS SEPARATED BY COMMAS'/) READ(ICC,4)ANS 4 FORMAT(7(A5,1X)) IF (ANS(1).EQ.'!') RETURN DO 5 J=1,7 IF (ANS(J).EQ.' ') GO TO 12 IF (ANS(J).NE.'HELP') GO TO 7 WRITE(IDLG,6) 6 FORMAT('0THE CHI SQUARE COMMAND OPERATES FROM RAW DATA NOT'/ 1 ' PRE CALCULATED TABLES, OPTIONS AVAILABLE ARE:'/ 2 ' "GROUP" - USED TO ESTABLISH GROUPINGS PRIOR TO COLPS'/ 3' "CONTG" - USED TO ELIMINATE CONTINGENCY TABLE FROM FINAL OUTPUT' 4/ ' "COLPS" - COLLAPSE CONTINGENCY TABLE (INCLUDES OMITTING)'/ 5' "FISHR" - FISHER''S EXACT TEST PROBABILITY (2X2 ONLY)'/ 6'0IF NO OPTIONS ARE DESIRED TYPE A RETURN') GO TO 1 7 IF (ANS(J).NE.'GROUP') GO TO 8 BREAK=1 GO TO 5 8 IF (ANS(J).NE.'CONTG') GO TO 9 CONTG=1 GO TO 5 9 IF (ANS(J).NE.'COLPS') GO TO 10 COLPS=1 GO TO 5 10 IF(ANS(J).NE.'FISHR') GO TO 19 FISHR=1 GO TO 5 19 WRITE(IDLG,11) ANS(J) 11 FORMAT(' OPTION "',A5,'" DOES NOT EXIST') GO TO 1 5 CONTINUE 12 IF(ICC.NE.2) WRITE(IDLG,13) 13 FORMAT(' WHICH VARIABLE OR VARIABLES? ',$) IRET=-98 CALL ALPHA(IV,2,NN,IRET,IHELP,IERR,NAMES,NV) IF (IRET.EQ.1) RETURN IF (IERR.EQ.1) GO TO 12 IF (IHELP.EQ.1) GO TO 12 IVB(1)=IV(1) IVB(2)=IV(2) IF (IV(1).LT.0) IVB(1)=1 IF (NN.LT.2) GO TO 18 IF (IV(2).LT.0) IVB(2)=1 IF ((IV(1).LT.0).AND.(IV(2).LT.0)) IVB(2)=2 GO TO 18 15 IF (IV(2).LT.0) GO TO 16 IF (IV(1).LT.0) GO TO 17 RETURN 16 IVB(2)=IVB(2)+1 IF (IVB(2).LE.NV) GO TO 18 IF(IV(1).GE.0) RETURN IVB(1)=IVB(1)+1 IVB(2)=IVB(1)+1 IF(IVB(2).LE.NV) GO TO 18 RETURN 17 IVB(1)=IVB(1)+1 IF (IVB(1).LE.NV) GO TO 18 RETURN 18 DO 2 I=1,2 2 LVL(I)=0 IF (NN.LT.2) GO TO 20 IF (IVB(1).EQ.IVB(2)) GO TO 15 C C HERE BEGINS ACTUAL CHI SQUARE C 20 IF (BREAK.NE.1) GO TO 35 DO 33 I=1,NN IF(ICC.NE.2) WRITE(IDLG,22) NAMES(IVB(I)) 22 FORMAT(' ENTER RANGES, 1 PER LINE FOR VARIABLE: ', A5/) DO 23 J=1,15 27 IF(ICC.NE.2) WRITE(IDLG,24) 24 FORMAT('+?',$) READ(ICC,25,END=21) ANS(1) 25 FORMAT(A5) IF (ANS(1).EQ.'!') RETURN IF (ANS(1).EQ.'STOP') GO TO 21 IF (ANS(1).EQ.' ') GO TO 21 IF(ANS(1).EQ.'AUTO') GO TO 32 IF (ANS(1).NE.'HELP') GO TO 28 WRITE(IDLG,26) 26 FORMAT(' ENTER EACH RANGE AFTER THE ? , SMALLER FIRST SEPARATED 1 BY A COMMA, ONE PER LINE'/ 2 ' EXAMPLE:'/' ?13,2'/' TO SIGNAL THE END OF RANGES TYPE 3 "STOP", ^Z, OR '/) GO TO 27 32 LVL(I)=0 GO TO 33 28 REREAD 29, (BRK(I,K,J),K=1,2) 29 FORMAT(2F) IF (BRK(I,1,J).LE.BRK(I,2,J)) GO TO 23 WRITE(IDLG,30) 30 FORMAT(' SMALLER FIRST'//) GO TO 27 23 BRK(I,3,J)=0 WRITE(IDLG,31) 31 FORMAT(' 15 RANGES MAXIMUM') J=16 21 LVL(I)=J-1 33 CONTINUE C C SET LEVELS FOR EACH SET OF VARIABLES C 35 DO 34 I=1,NC ISUB(I)=I 34 LEVEL(I)=0 NK=NC DO 36 I=1,NN KONS=1 IF(I.EQ.2) KONS=IJUST IF (LVL(I).EQ.0) GO TO 39 DO 37 J=1,NK DO 38 K=1,LVL(I) IF (DATA(ISUB(J),IVB(I)).LT.BRK(I,1,K)) GO TO 38 IF (DATA(ISUB(J),IVB(I)).GT.BRK(I,2,K)) GO TO 38 BRK(I,3,K)=1 LEVEL(J)=K*KONS+LEVEL(J) GO TO 37 38 CONTINUE ISUB(J)=0 37 CONTINUE DO 130 J=1,LVL(I) IF(BRK(I,3,J).NE.0) GO TO 130 WRITE(IDLG,131) BRK(I,1,J),BRK(I,2,J) 131 FORMAT(' RANGE ',G10.4,',',G10.4,' HAS NO VALUES - LEVEL DROPPED') DO 132 K=1,NK IF(LEVEL(K).GT.J) LEVEL(K)=LEVEL(K)-1*KONS 132 CONTINUE 130 CONTINUE GO TO 41 39 LL=I+2 CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,LL) L=1 VALUE=DATA(ISUB(1),IVB(I)) DO 40 J=1,NK IF (DATA(ISUB(J),IVB(I)).EQ.VALUE) GO TO 40 VALUE=DATA(ISUB(J),IVB(I)) L=L+1 40 LEVEL(J)=L*KONS+LEVEL(J) LVL(I)=L 41 J=0 K=1 42 IF (ISUB(K).EQ.0) GO TO 43 J=J+1 ISUB(J)=ISUB(K) LEVEL(J)=LEVEL(K) 43 K=K+1 IF (K.LE.NK) GO TO 42 NK=J 36 CONTINUE IF(NK.LT.1) RETURN C C LEVELS DETERMINED NOW C IF (COLPS.NE.1) GO TO 100 C C COLLAPSING PORTION C CALL CHITAB(NK,LEVEL,ISUB,MV,MC,DATA,NAMES,IVB,1,NN,LINES) IF(ICC.NE.2) WRITE(IDLG,50) 50 FORMAT('0COLLAPSING PORTION INSERT 1 INSTRUCTION PER LINE AFTER 1 THE ?') 51 IF(ICC.NE.2) WRITE(IDLG,52) 52 FORMAT(' ? ',$) READ(ICC,53,END=95) STRG 53 FORMAT(50A1) IF (STRG(1).EQ.' ') GO TO 95 IF (STRG(1).EQ.'!') RETURN IF ((STRG(1).EQ.'S').AND.(STRG(2).EQ.'T').AND.(STRG(3).EQ.'O') 1.AND.(STRG(4).EQ.'P')) GO TO 95 IF ((STRG(1).NE.'H').OR.(STRG(2).NE.'E').OR.(STRG(3).NE.'L') 1.OR.(STRG(4).NE.'P')) GO TO 55 WRITE(IDLG,54) 54 FORMAT('0TWO OPERATIONS ARE AVAILABLE TO COLLAPSE THE CHI-SQUARE:'/ 1'"C"-COMBINE AND "D"-COLLAPSE. THEY MUST BE THE FIRST CHARACTER 2 OF'/' THE LINE. THIS MUST BE FOLLOWED BY THE VARIABLE NAME 3 OR NUMBER'/' ENCLOSED IN PARENTHESIS. THIS IS FOLLOWED BY 4 A STRING OF NUMBERS'/' REFERENCING LEVELS OF THE VARIABLES. 5 RANGES MAY BE SPECIFIED'/' BY SEPARATEING THE EXTREMES OF', 5' THE RANGE WITH A "-".'/' A , ^Z, OR "STOP"', 6' WILL EXIT FROM THE'/' COLLAPSING PORTION AND CONTINUE THE 7 CHI SQUARE'/' EXAMPLE:'/' C(AGE)1-5'/' COMBINE THE LEVELS 8 1 THRU 5 FOR THE VARIABLE AGE.'/' D(3)1,9,3'/' DELETE 9 LEVELS 1,9,AND 3 FOR VARIABLE 3') GO TO 51 55 IF ((STRG(1).EQ.'C').OR.(STRG(1).EQ.'D')) GO TO 57 WRITE(IDLG,56) 56 FORMAT(' INSTRUCTION MUST BE "C" OR "D"') GO TO 51 57 IF (STRG(2).EQ.'(') GO TO 60 58 WRITE(IDLG,59) 59 FORMAT(' VARIABLE NOT SPECIFIED, OR NOT IN PARENTHESIS') GO TO 51 C C DETERMINE VARIABLE TO BE USED IN COLLAPSING C 60 DO 61 I=1,5 61 B(I)=' ' NUM=0 IF ((STRG(3).LE.'9').AND.(STRG(3).GE.'0')) NUM=1 K=1 J=3 62 IF (STRG(J).EQ.')') GO TO 65 IF (NUM.EQ.0) GO TO 63 IF ((STRG(J).LE.'9').AND.(STRG(J).GE.'0')) GO TO 63 GO TO 58 63 IF (K.GT.5) GO TO 64 B(K)=STRG(J) K=K+1 64 J=J+1 IF (J.LE.50) GO TO 62 65 IF (NUM.EQ.0) GO TO 165 IF (B(5).NE.' ') GO TO 165 DO 166 I=4,1,-1 166 B(I+1)=B(I) B(1)=' ' GO TO 65 165 ENCODE(5,66,NNS) B 66 FORMAT(5A1) IF (NUM.EQ.0) GO TO 70 DECODE(5,67,NNS) NS 67 FORMAT(I5) 72 IF ((IVB(1).EQ.NS).OR.(IVB(2).EQ.NS)) GO TO 75 68 WRITE(IDLG,69) 69 FORMAT(' VARIABLE SPECIFIED IS NOT BEING ANALYZED') GO TO 51 70 DO 71 NS=1,NV IF (NAMES(NS).EQ.NNS) GO TO 72 71 CONTINUE WRITE(IDLG,73) NNS 73 FORMAT(' VARIABLE "',A5,'" DOES NOT EXIST') GO TO 51 C C DETERMINE THE LEVELS TO BE ACTED UPON C 75 K=1 IF (NS.EQ.IVB(2)) K=2 THRU=0 M=1 175 J=J+1 L=1 DO 76 I=1,5 76 B(I)=' ' 78 IF (STRG(J).EQ.' ') GO TO 79 IF (STRG(J).EQ.',') GO TO 79 IF (STRG(J).EQ.'-') GO TO 79 IF ((STRG(J).LT.'0').OR.(STRG(J).GT.'9')) GO TO 181 IF (L.GT.5) GO TO 77 B(L)=STRG(J) L=L+1 77 J=J+1 IF (J.LE.50) GO TO 78 79 IF (L.EQ.1) GO TO 83 80 IF (B(5).NE.' ') GO TO 82 DO 81 I=4,1,-1 81 B(I+1)=B(I) B(1)=' ' GO TO 80 181 WRITE(IDLG,182) 182 FORMAT(' LEVELS MAY ONLY BE SPECIFIED BY LEVEL NUMBERS') GO TO 51 82 ENCODE(5,66,NNS) B DECODE(5,67,NNS)NS IF ((NS.LE.0).OR.(NS.GT.LVL(K))) GO TO 183 LTBM(M)=NS IF (THRU.EQ.1) LTBM(M)=-NS THRU=0 M=M+1 IF (STRG(J).EQ.'-') THRU=1 IF (STRG(J).EQ.' ') GO TO 83 GO TO 175 183 WRITE(IDLG,184) LVL(K) 184 FORMAT(' LEVEL SPECIFIED DOES NOT EXIST, MAXIMUM LEVEL IS',I4) GO TO 51 83 IF (M.GT.1) GO TO 85 WRITE(IDLG,84) 84 FORMAT(' NO LEVELS TO BE ACTED UPON') GO TO 51 85 M=M-1 KMIN=LTBM(1) DO 86 I=2,M L=LTBM(I) IF (L.LT.0) L=-L IF (KMIN.LE.L) GO TO 87 KMIN=L 87 IF (LTBM(I).GT.0) GO TO 86 L=-LTBM(I) IF (LTBM(I-1).GT.0) GO TO 89 WRITE(IDLG,88) 88 FORMAT(' TWO - IN A ROW') GO TO 51 89 IF (LTBM(I-1).LT.L) GO TO 86 LTBM(I)=-LTBM(I-1) LTBM(I-1)=L 86 CONTINUE C C ACTUALLY GO THRU AND RECODE THOSE VARIABLES USED IN COLAPSE C PORTION C J=1 190 L=LEVEL(J).AND.IRMAP IF(K.EQ.2) L=(LEVEL(J).AND.ILMAP)/IJUST DO 91 I=1,M IF (I.GE.M) GO TO 92 IF (LTBM(I+1).GT.0) GO TO 92 IF ((L.LT.LTBM(I)).OR.(L.GT.-LTBM(I+1))) GO TO 91 GO TO 93 92 IF (L.NE.LTBM(I)) GO TO 91 93 IF(STRG(1).EQ.'D') GO TO 193 IF(K.EQ.2) LEVEL(J)=(IRMAP.AND.LEVEL(J))+KMIN*IJUST IF(K.EQ.1) LEVEL(J)=(ILMAP.AND.LEVEL(J))+KMIN GO TO 90 193 IF(J.EQ.NK) GO TO 194 DO 94 II=J,NK-1 ISUB(II)=ISUB(II+1) 94 LEVEL(II)=LEVEL(II+1) 194 NK=NK-1 IF(J.GT.NK) GO TO 51 GO TO 190 91 CONTINUE 90 J=J+1 IF(J.LE.NK) GO TO 190 GO TO 51 C C RECODE LEVELS AFTER COLLAPSING C 95 DO 96 I=1,NN L=I CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,L) L=1 KVAL=LEVEL(1).AND.IRMAP IF(I.EQ.2) KVAL=(LEVEL(1).AND.ILMAP)/IJUST DO 97 J=1,NK K=LEVEL(J).AND.IRMAP IF(I.EQ.2) K=(LEVEL(J).AND.ILMAP)/IJUST IF (K.EQ.KVAL) GO TO 98 L=L+1 KVAL=K 98 IF(I.EQ.2) LEVEL(J)=(IRMAP.AND.LEVEL(J)) +L*IJUST IF(I.EQ.1) LEVEL(J)=(ILMAP.AND.LEVEL(J))+L 97 CONTINUE 96 LVL(I)=L IF(NK.LT.1) RETURN C C BEGIN OUTPUT PHASE C 100 IF(IOUT.NE.21)WRITE(IOUT,5566) (HEDR(I),I=1,NSZ) 5566 FORMAT('1',70A1) IF(IOUT.EQ.21) CALL PRNTHD WRITE(IOUT,103) 103 FORMAT('0',20X,'***** CHI SQUARE *****') LINES=4 IF (CONTG.EQ.1) GO TO 101 CALL CHITAB(NK,LEVEL,ISUB,MV,MC,DATA,NAMES,IVB,2,NN,LINES) 101 IF (NN.GT.1) GO TO 105 110 A=LVL(1) IF(A.EQ.0) GO TO 106 EXP=NK/A CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,5) CHI=0 KVAL=LEVEL(1) O=1 DO 102 J=2,NK IF (KVAL.EQ.LEVEL(J)) GO TO 102 CHI=CHI+(O-EXP)**2/EXP O=0 KVAL=LEVEL(J) 102 O=O+1 CHI=CHI+(O-EXP)**2/EXP NF=LVL(1)-1 IF (NF.LT.1) GO TO 106 IF(IOUT.NE.21) GO TO 111 LINES=LINES+2 IF(LINES.LE.(LINPP-1)) GO TO 111 CALL PRNTHD LINES=4 111 WRITE(IOUT,104) CHI,NF 104 FORMAT('0CHI SQUARE =',G,' WITH',I4,' DEGREES OF FREEDOM') GO TO 123 C C CHECK TO SEE IF IT MAY BE ANALY. AS A ONE SAMPLE CHI SQ. C 106 IF(IOUT.NE.21) GO TO 112 LINES=LINES+2 IF(LINES.LE.LINPP) GO TO 112 CALL PRNTHD LINES=4 112 WRITE(IOUT, 107) 107 FORMAT('0CONTINGENCY TABLE DEGENERATED TO 1 CELL - NO CHI 1 SQUARE PERFORMED') GO TO 15 108 IF (LVL(2).GT.LVL(1)) LVL(1)=LVL(2) LVL(2)=0 WRITE(IDLG,109) 109 FORMAT(' CHI SQUARE ANALYSIS PERFORMED FOR 1 SAMPLE CASE') GO TO 110 105 IF(LVL(1).LT.2) GO TO 108 IF(LVL(2).LT.2) GO TO 108 C C TWO SAMPLE CHI SQUARE C NF=(LVL(1)-1)*(LVL(2)-1) CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,5) C C BREAK IT INTO TWO GROUPS FOR CALCULATION OF CHI SQUARE C DO 115 I=1,NK ISUB(I)=(LEVEL(I).AND.ILMAP)/IJUST 115 LEVEL(I)=LEVEL(I).AND.IRMAP C C CALCULATION FOR CHI SQUARE C CHI=0 K=1 DO 116 I=1,LVL(2) RO=0 117 IF (ISUB(K).NE.I) GO TO 120 RO=RO+1 K=K+1 IF (K.LE.NK) GO TO 117 120 DO 118 J=1,LVL(1) CO=0 OBS=0 DO 119 L=1,NK IF (LEVEL(L).NE.J) GO TO 119 CO=CO+1 IF (ISUB(L).EQ.I) OBS=OBS+1 119 CONTINUE EXP=(CO*RO)/NK 118 CHI=CHI+(OBS-EXP)**2/EXP 116 CONTINUE IF(IOUT.NE.21) GO TO 113 LINES=LINES+2 IF(LINES.LE.(LINPP-1)) GO TO 113 CALL PRNTHD LINES=4 113 WRITE(IOUT,104) CHI,NF IF (NF.NE.1) GO TO 123 C C YATES CORRECTION FACTOR AUTOMATICALLY PERFORMED IF DEGREES OF C FREEDOM IS EQUAL TO ONE C DO 125 I=1,2 DO 125 J=1,2 125 MCC(I,J)=0 DO 121 I=1,NK 121 MCC(ISUB(I),LEVEL(I))=MCC(ISUB(I),LEVEL(I))+1. CHI=MCC(1,1)*MCC(2,2)-MCC(1,2)*MCC(2,1) IF (CHI.LT.0) CHI=-CHI CHI=(CHI-NK/2.)**2 CHI=(CHI*NK)/((MCC(1,1)+MCC(1,2))*(MCC(1,1)+MCC(2,1))*(MCC(2,2) 1+MCC(1,2))*(MCC(2,2)+MCC(2,1))) IF(IOUT.NE.21) GO TO 128 LINES=LINES+2 IF(LINES.LE.(LINPP-1)) GO TO 128 CALL PRNTHD LINES=4 128 WRITE(IOUT,122) CHI 122 FORMAT('0CORRECTED CHI SQUARE =',G) 123 CHI=CHI/NF PROB=FISHER(NF,1000,CHI) IF(IOUT.NE.21) GO TO 129 LINES=LINES+1 IF(LINES.LE.LINPP) GO TO 129 CALL PRNTHD LINES=3 129 WRITE(IOUT,124) PROB 124 FORMAT(' HAVING A PROBABILITY OF ',F5.2) IF((NF.NE.1).OR.(FISHR.EQ.0).OR.(LVL(2).EQ.0)) GO TO 137 C C SECTION TO CALCULATE FISHERS EXACT PROBABILITY C DO 138 I=1,2 MCC(I,3)=0 138 MCC(3,I)=0 DO 141 I=1,2 MCC(I,3)=MCC(I,1)+MCC(I,2) 141 MCC(3,I)=MCC(1,I)+MCC(2,I) MCC(3,3)=MCC(1,3)+MCC(2,3) CALL CONP(MCC,3,3,PT,PS,PC) IF(IOUT.NE.21) GO TO 142 LINES=LINES+4 IF(LINES.LE.LINPP) GO TO 142 CALL PRNTHD LINES=6 142 WRITE(IOUT,127) PT,PS 127 FORMAT('0FISHER''S EXACT PROBABILITY FOR OBTAINING THE ', 1'GIVEN TABLE =',F8.5/23X,'FOR OBTAINING A TABLE AS ', 2'PROBABLE, OR'/23X,'LESS PROBABLE THAN THE GIVEN', 3' TABLE =',F8.5) C C 137 IF(NN.EQ.1) GO TO 15 CONCOF=SQRT(CHI/(NK+CHI)) IF(IOUT.NE.21) GO TO 143 LINES=LINES+2 IF(LINES.LE.LINPP) GO TO 143 CALL PRNTHD LINES=4 143 WRITE(IOUT,126) CONCOF 126 FORMAT('0CONTINGENCY COEFFICIENT =',G) GO TO 15 END C *** STAT PACK *** C SUBROUTINE CALLED BY CHI TO WRITE OUT CONTINGENCY TABLE C CALLING SEQUENCE: CALL CHITAB(NL,LEVEL,ISUB,MV,MC,DATA,NAMES,IX,ITM,NUMV,LINES) C WHERE NL - IS THE NUMBER OF OBS. IN CONTG TABLE C LEVEL - IS A VECTOR DIMENSIONED AT LEAST NL C ISUB - IS A VECTOR DIMENSIONED AT LEAST NL C MV - MAX. NUMBER OF VARIABLES C MC - MAX. NUMBER OF OBSERVATIONS C DATA - MATRIX FOR DATA, DIMENSIONED FOR MAXIMUM C NAMES - A VECTOR CONTAINING VARIABLE NAMES C ITM - SWITCH TO SHOW HOW THINGS SHOULD BE OUTPUT C NUMV - ONE OR TWO WAY CHISQ? C LINES - NUMBER OF LINES PRESENTLY USED ON THE LPT PAGE C C SIMPLY TAKES DATA REFERENCED BY ISUB, LEVELS DETERMINED C BY LEVEL AND CREATES A CONTINGENCY TABLE C SUBROUTINE CHITAB (NL,LEVEL,ISUB,MV,MC,DATA,NAMES,IX,ITM,NUMV, 1LINES) COMMON/DEV/ICC,IDATA,IOUT,IDLG,IDSK COMMON /PRNT/ LINPP,ICOPS,RUNPRG DIMENSION LEVEL(1), ISUB(1), DATA(MC,MV), NAMES(1), XOUT(20), 1IX(2), ICT(20), ITOT(20),IHDR(20) ILMAP="007777700000 IRMAP="000000077777 IJUST="000000100000 ILL=10 IF(ITM.EQ.1) GO TO 8 IF(IOUT.EQ.21) ILL=20 8 DOTS='.....' ITAB0=IOUT IF(ITM.EQ.1) ITAB0=IDLG DO 10 KKK=1,NUMV KKKS=KKK CALL CHISOR(NL,LEVEL,ISUB,IX,DATA,MV,MC,KKKS) IF(ITAB0.NE.21) GO TO 40 LINES=LINES+2 IF(LINES.LE.(LINPP-10)) GO TO 40 CALL PRNTHD LINES=4 40 IF((KKK.EQ.2).AND.(NUMV.NE.1)) WRITE(ITAB0,1) NAMES(IX(KKK)) 1 FORMAT ('0LEVELS FOR VERTICAL VARIABLE: ',A5) IF((KKK.EQ.1).AND.(NUMV.NE.1)) WRITE(ITAB0,7) NAMES(IX(KKK)) 7 FORMAT ('0LEVELS FOR HORIZONTAL VARIABLE: ',A5) IF(NUMV.EQ.1) WRITE(ITAB0,30)NAMES(IX(1)) 30 FORMAT('0LEVELS FOR VARIABLE: ',A5) WRITE(ITAB0,2) 2 FORMAT ('0LEVEL',2X,'VALUES COMPRISING LEVEL') LEVL=1 VALUE=DATA(ISUB(1),IX(KKK)) IFRQ=1 K=1 XOUT(K)=VALUE K=K+1 ISW=0 DO 3 I=2,NL LH=LEVEL(I).AND.IRMAP IF(KKK.EQ.2) LH=(LEVEL(I).AND.ILMAP)/IJUST IF(LEVL.NE.LH) GO TO 4 IF(VALUE.EQ.DATA(ISUB(I),IX(KKK))) GO TO 3 XOUT(K)=',' K=K+1 IF(K.GT.ILL) GO TO 4 VALUE=DATA(ISUB(I),IX(KKK)) XOUT(K)=VALUE K=K+1 GO TO 3 4 IF(ITAB0.NE.21) GO TO 41 LINES=LINES+1 IF(ISW.EQ.0) LINES=LINES+1 IF(LINES.LE.(LINPP-1)) GO TO 41 CALL PRNTHD LINES=3 IF(ISW.EQ.0) LINES=LINES+1 41 IF(ISW.EQ.0) WRITE(ITAB0,5) LEVL,(XOUT(L),L=1,K-1) 5 FORMAT ('0',I3,4X,10(G10.4,A2)) IF(ISW.EQ.1) WRITE (ITAB0,6)(XOUT(L),L=1,K-1) 6 FORMAT (8X,10(G10.4,A2)) ISW=1 K=1 VALUE=DATA(ISUB(I),IX(KKK)) XOUT(K)=VALUE K=K+1 IF(LEVL.EQ.LH) GO TO 3 IF(NUMV.EQ.1) WRITE(ITAB0,31) IFRQ 31 FORMAT(9X,'FREQUENCY =',I4) IFRQ=0 ISW=0 LEVL=LEVL+1 GO TO 3 3 IFRQ=IFRQ+1 IF(ITAB0.NE.21) GO TO 43 LINES=LINES+1 IF(ISW.EQ.0) LINES=LINES+1 IF(LINES.LE.(LINPP-1)) GO TO 43 CALL PRNTHD LINES=3 IF(ISW.EQ.0) LINES=LINES+1 43 IF(ISW.EQ.0) WRITE(ITAB0,5) LEVL,(XOUT(L),L=1,K-1) IF(ISW.EQ.1) WRITE(ITAB0,6)(XOUT(L),L=1,K-1) IF(NUMV.EQ.1) WRITE(ITAB0,31)IFRQ IF(KKK.EQ.1) MAX=LEVL+1 10 CONTINUE IF(NUMV.EQ.1) RETURN LL=MAX/ILL IF((MAX-ILL*LL).NE.0) LL=LL+1 DO 11 I=1,LL IADD=(I-1)*ILL IEND=ILL IF((IEND+IADD).GT.MAX) IEND=MAX-IADD DO 12 J=1,IEND L=IADD+J ENCODE(5,20,IHDR(J)) L 20 FORMAT (I5) ITOT(J)=0 12 CONTINUE IF((IADD+IEND).EQ.MAX) IHDR(IEND)='TOTAL' IF(ITAB0.NE.21) GO TO 45 LINES=LINES+4 IF(LINES.LE.(LINPP-10)) GO TO 45 CALL PRNTHD LINES=6 45 WRITE (ITAB0,13) NAMES(IX(1)), (IHDR(J),J=1,IEND) WRITE (ITAB0,14) NAMES(IX(2)),(DOTS,DOTS,J=1,IEND) 13 FORMAT ('0',30X,A5/' LEVEL',2X,20(1X,A5)) 14 FORMAT (1X,A5,1X,'.',40A3) LEVL=1 DO 15 J=1,IEND 15 ICT(J)=0 K=1 23 LH=(LEVEL(K).AND.ILMAP)/IJUST IRH=LEVEL(K).AND.IRMAP IF(LEVL.NE.LH) GO TO 17 IF((IEND+IADD).EQ.MAX) ICT(IEND)=ICT(IEND)+1 IF((IRH.LT.IADD+1).OR.(IRH.GT.IADD+IEND)) GO TO 16 ICT(IRH-IADD)=ICT(IRH-IADD)+1 GO TO 16 17 IF(ITAB0.NE.21) GO TO 46 LINES=LINES+1 IF(LINES.LE.LINPP) GO TO 46 CALL PRNTHD WRITE (ITAB0,13) NAMES(IX(1)), (IHDR(J),J=1,IEND) WRITE (ITAB0,14) NAMES(IX(2)),(DOTS,DOTS,J=1,IEND) LINES=7 46 WRITE (ITAB0,18) LEVL,(ICT(J),J=1,IEND) 18 FORMAT (1X,I3,3X,'.',20(1X,I5)) DO 19 J=1,IEND ITOT(J)=ITOT(J)+ICT(J) 19 ICT(J)=0 LEVL=LEVL+1 GO TO 23 16 K=K+1 IF(K.LE.NL) GO TO 23 DO 22 J=1,IEND 22 ITOT(J)=ITOT(J)+ICT(J) IF(ITAB0.NE.21) GO TO 47 LINES=LINES+2 IF(LINES.LE.LINPP) GO TO 47 CALL PRNTHD WRITE (ITAB0,13) NAMES(IX(1)), (IHDR(J),J=1,IEND) WRITE (ITAB0,14) NAMES(IX(2)),(DOTS,DOTS,J=1,IEND) LINES=8 47 WRITE(ITAB0,18) LEVL,(ICT(J),J=1,IEND) WRITE (ITAB0,21)(ITOT(J),J=1,IEND) 21 FORMAT (1X,'TOTAL',1X,'.',20(1X,I5)) 11 CONTINUE RETURN END C *** STAT PACK *** C SUBROUTINE PART OF CHI USED TO SORT C CALLING SEQUENCE: CALL CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,ITYP) C WHRER NK - IS THE NUMBER OF OBSERVATIONS TO BE SORTED C LEVEL IS A VECTOR DIMENSIONED AT LEAST NK C ISUB IS A VECTOR DIMENSIONED AT LEAST NK C IBV - IS A VECTOR CONTAINING THE VARIABLE NUMBERS TO BE USED C DATA - IS MATRIX CONTAINING DATA C MV - IS MAXIMUM NUMBER OF VARIABLES C MC - IS MAXIMUM NUMBER OF OBSERVATIONS C ITYP - TELLS IN WHICH WAY THEY ARE TO BE SORTED C C SORTS NUMBERS BY SEVERAL DIFFERENT SEQUENCES C SUBROUTINE CHISOR(NK,LEVEL,ISUB,IVB,DATA,MV,MC,ITYP) DIMENSION DATA(MC,MV),LEVEL(1),ISUB(1),IVB(2) DIMENSION IU(16),IL(16) ILMAP="007777700000 IRMAP="000000077777 IBOTH="007777777777 IF((ITYP.NE.3).AND.(ITYP.NE.4)) GO TO 1 IFIN=0 IS=ITYP-2 GO TO 2 1 IF(ITYP.EQ.1) IFIN=IRMAP IF(ITYP.EQ.2) IFIN=ILMAP IF(ITYP.EQ.5) IFIN=IBOTH IS=ITYP IF(IS.GT.2) IS=1 2 IS=IVB(IS) M=1 II=1 J=NK 71 IF(II.GE.J) GO TO 78 72 K=II IJ=(J+II)/2 T=DATA(ISUB(IJ),IS) IT=IFIN.AND.LEVEL(IJ) KT=LEVEL(II).AND.IFIN IF(KT.LT.IT) GO TO 73 IF((KT.EQ.IT).AND.(DATA(ISUB(II),IS).LE.T)) GO TO 73 ISAV=LEVEL(IJ) LEVEL(IJ)=LEVEL(II) LEVEL(II)=ISAV ISAV=ISUB(IJ) ISUB(IJ)=ISUB(II) ISUB(II)=ISAV IT=KT T=DATA(ISUB(IJ),IS) 73 LL=J KT=LEVEL(J).AND.IFIN IF(KT.GT.IT) GO TO 75 IF((KT.EQ.IT).AND.(DATA(ISUB(J),IS).GE.T)) GO TO 75 ISAV=LEVEL(IJ) LEVEL(IJ)=LEVEL(J) LEVEL(J)=ISAV ISAV=ISUB(IJ) ISUB(IJ)=ISUB(J) ISUB(J)=ISAV IT=KT T=DATA(ISUB(IJ),IS) KT=LEVEL(II).AND.IFIN IF(KT.LT.IT) GO TO 75 IF((KT.EQ.IT).AND.(DATA(ISUB(II),IS).LE.T)) GO TO 75 ISAV=LEVEL(IJ) LEVEL(IJ)=LEVEL(II) LEVEL(II)=ISAV ISAV=ISUB(IJ) ISUB(IJ)=ISUB(II) ISUB(II)=ISAV IT=KT T=DATA(ISUB(IJ),IS) GO TO 75 74 ISAV=LEVEL(LL) LEVEL(LL)=LEVEL(K) LEVEL(K)=ISAV ISAV=ISUB(LL) ISUB(LL)=ISUB(K) ISUB(K)=ISAV 75 LL=LL-1 KT=LEVEL(LL).AND.IFIN IF(KT.GT.IT) GO TO 75 IF((KT.EQ.IT).AND.(DATA(ISUB(LL),IS).GT.T)) GO TO 75 76 K=K+1 KT=LEVEL(K).AND.IFIN IF(KT.LT.IT) GO TO 76 IF((KT.EQ.IT).AND.(DATA(ISUB(K),IS).LT.T)) GO TO 76 IF(K.LE.LL) GO TO 74 IF((LL-II).LE.(J-K)) GO TO 77 IL(M)=II IU(M)=LL II=K M=M+1 GO TO 79 77 IL(M)=K IU(M)=J J=LL M=M+1 GO TO 79 78 M=M-1 IF(M.EQ.0) GO TO 90 II=IL(M) J=IU(M) 79 IF((J-II).GE.11) GO TO 72 IF(II.EQ.1) GO TO 71 II=II-1 80 II=II+1 IF(II.EQ.J) GO TO 78 T=DATA(ISUB(II+1),IS) IT=LEVEL(II+1).AND.IFIN KT=LEVEL(II).AND.IFIN IF(KT.LT.IT) GO TO 80 IF((KT.EQ.IT).AND.(DATA(ISUB(II),IS).LE.T)) GO TO 80 K=II ISAV=LEVEL(II+1) ISAV2=ISUB(II+1) 81 LEVEL(K+1)=LEVEL(K) ISUB(K+1)=ISUB(K) K=K-1 KT=LEVEL(K).AND.IFIN IF(IT.LT.KT) GO TO 81 IF((IT.EQ.KT).AND.(T.LT.DATA(ISUB(K),IS))) GO TO 81 LEVEL(K+1)=ISAV ISUB(K+1)=ISAV2 GO TO 80 90 RETURN END C *** STAT PACK *** C SUBROUTINE IS PART OF STAT PACK CHISQUARE USED TO CALCULATE C THE EXACT PROBABILITY OF A 2X2 TABLE OR THE PROBABILITY OF C HAVING THAT TABLE OR A TABLE LESS PROBABLE THAN IT. C ROUTINE FROM COMMUNICATIONS OF ACM NOVEMBER 1972 C WHERE MATRIX - IS A 3X3 TABLE FOR WHICH THE PROB IS FOUND C NR - IS THE NUMBER OF ROWS (HERE ALWAYS 2) C NC - IS THE NUMBER OF COLUMNS (HERE ALWAYS 2) C PT - PROBABILITY OF GIVEN TABLE C PS - PROBABILITY OF TABLE AS PROBABLE AS GIVEN TABLE C PC - NOT USED IN STP, BUT IT IS THE PROBABILITY OF C OBTAINING SOME OF THE TABLES POSSIBLE WITHIN C THE CONSTRAINTS OF THE MARGINAL TOTALS C C SUBROUTINE CONP(MATRIX,NR,NC,PT,PS,PC) DIMENSION MATRIX(NR,NC) COMMON/IFLAG1/IFLAG INTEGER R,C,TEMP IFLAG=0 R=NR-1 C=NC-1 QXLOG=-FACLOG(MATRIX(NR,NC)) DO 10 I=1,R 10 QXLOG=QXLOG+FACLOG(MATRIX(I,NC)) DO 20 J=1,C 20 QXLOG=QXLOG+FACLOG(MATRIX(NR,J)) RXLOG=0.0 DO 50 I=1,R DO 50 J=1,C 50 RXLOG=RXLOG+FACLOG(MATRIX(I,J)) PT=10.0**(QXLOG-RXLOG) PS=0 PC=0 DO 100 I=2,R DO 100 J=2,C 100 MATRIX(I,J)=MIN0(MATRIX(I,NC),MATRIX(NR,J)) GO TO 300 200 DO 220 I=2,R DO 220 J=2,C MATRIX(I,J)=MATRIX(I,J)-1 IF(MATRIX(I,J).GE.0) GO TO 300 220 MATRIX(I,J)=MIN0(MATRIX(I,NC),MATRIX(NR,J)) RETURN 300 DO 320 I=2,R TEMP=MATRIX(I,NC) DO 310 J=2,C 310 TEMP=TEMP-MATRIX(I,J) IF(TEMP.LT.0) GO TO 200 320 MATRIX(I,1)=TEMP DO 340 J=1,C TEMP=MATRIX(NR,J) DO 330 I=2,R 330 TEMP=TEMP-MATRIX(I,J) IF(TEMP.LT.0) GO TO 200 340 MATRIX(1,J)=TEMP RXLOG=0.0 DO 350 I=1,R DO 350 J=1,C 350 RXLOG=RXLOG+FACLOG(MATRIX(I,J)) PX=10.0**(QXLOG-RXLOG) PC=PC+PX IF((PT/PX).GT.0.9999) PS=PS+PX GO TO 200 END C *** STAT PACK *** C FUNCTION IS PART OF STP TAKEN FOR COMMUNICATIONS OF ACM C NOVEMBER 1972. USED TOR FINDING LOG BASE 10 OF N FACTORIAL. C USES STIRLINGS APPROX. IF N.GT.100 C C FIRST TIME THRU IT CREATES A TABLE IT USES. IFLAG IS THE C INDICATOR USED TO ESTABLISH RATHER A TABLE NEED BE CREATED C FUNCTION FACLOG(N) DIMENSION TABLE (101) COMMON /IFLAG1/IFLAG TPILOG=.3990899342 ELOG=.4342944819 IF(N.GT.100) GO TO 50 IF(IFLAG.EQ.0) GOTO 100 10 FACLOG=TABLE(N+1) RETURN 50 X=FLOAT(N) FACLOG=(X+.5)*ALOG10(X)-X*ELOG+TPILOG+ELOG/(12.0*X)-ELOG 1/(360.0*X**3) RETURN C C CREATE TABLE TO BE USED FOR REST OF TIME C 100 TABLE(1)=0.0 DO 120 I=2,101 X=FLOAT(I-1) 120 TABLE(I)=TABLE(I-1)+ALOG10(X) IFLAG=1 GOTO 10 END