C                                          *** STAT PACK ***
C     SUBROUTINE FOR 2 WAY ANALYSIS OF VARIANCE.
C     CALLING SEQUENCE: CALL ANOV2(NV,NC,MV,MC,DATA,VMN,STD,NAMES)
C     WHERE NV - IS THE NUMBER OF COLUMNS ACTUALLY FILLED (VARIABLES)
C           NC - IS THE NUMBER OF ROWS ACTUALLY FILLED (CASES)
C           MV - IS THE MAXIMUM NUMBER OF COLUMNS.
C           MC - IS THE MAXIMUM NUMBER OF ROWS.
C           DATA - DATA MATRIX, DIMENSIONED FOR MAXIMUM.
C           VMN - IS A VECTOR CONTAINING THE MEANS.
C           STD - IS A VECTOR CONTAINING THE STANDARD DEVIATIONS.
C           NAMES - IS A VECTOR CONTAINING VARIABLE NAMES
C
C     SUBROUTINE ALLOWS FOR EACH CELL TO BE AN INDIVIDUAL VARIABLE
C     OR FOR A BREAKDOWN ON THE VALUE OF ANOTHER VARIABLE TO DETERMINE
C     THE CELLS.  THIS PARTICULAR ANOV IS USED SINCE IT DOES ALLOW
C     EMPTY CELLS OR AN UNBALANCED DESIGN.  ORIGINALLY WRITTEN
C     AT W. M. U. BY EVA GAINES, UNDER DIRECTION OF MICHEAL STOLINE,
C     IT WAS MODIFIED DURING THE SUMMER OF 71 BY SAM ANEMA,
C     THE VERSION NOW IN STP IS A FURTHUR MODIFICATION OF THAT
C     PROGRAM.  ORIGINAL VERSION IS BASED ON BOOK BY BANCROFT
C     "INTERMEDIATE STATISTICAL METHODS".
C
      SUBROUTINE ANOV2(NV,NC,MV,MC,DATA,VMN,STD,NAMES)
      DIMENSION DATA(MC,MV), VMN(1), STD(1), OPT(5), OPTD(5)
      COMMON /DEV/ ICC, IDATA, IOUT, IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON/EXTRA/HEDR(70),NSZ
      DIMENSION YS(20,20),Y1(20),Y2(20),S1(20),S2(20),B1(20)
      DIMENSION T(20,20),C(20),AM(20),YT(20,20),IBK(2),BRNG(2,20,2)
      DIMENSION NGRPS(2),IVA(120)
      DIMENSION ICELLS(40),NAMES(1)
      EQUIVALENCE (NGRPS(1),NGRPA),(NGRPS(2),NGRPB)
      EQUIVALENCE (IBK(1),IBK1),(IBK(2),IBK2)
      EQUIVALENCE (IVA,ICELLS),(IVA(41),BRNG)
      ISQ=7
      IF(IOUT.EQ.21) ISQ=15
      ALL=0
5     IF(ICC.NE.2) WRITE (IDLG,1)
      OPTD(1)=0
      OPTD(2)=0
      OPTD(3)=0
      OPTD(4)=0
      OPTD(5)=0
      NGRPS(1)=20
      NGRPS(2)=20
1     FORMAT ('0LIST OPTIONS SEPARATED BY COMMAS'/)
      READ (ICC,3) OPT
2     FORMAT ('0THE 2 WAY ANALYSIS OF VARIANCE ASSUMES
     1 CELLS TO '/' BE MADE UP OF INDIVIDUAL VARIABLES
     2.'/' IT IS POSSIBLE HOWEVER TO BREAK
     3 1 VARIABLE INTO'/' CELLS, BY MEANS OF
     4 TWO OTHER VARIABLES.'/' IN THIS CASE, RANGES
     5 FOR BREAK-DOWNS ARE ASSUMED'/' TO BE SUPPLIED
     6 BY THE USER.  WHEN THE PROGRAM'/' CALLS FOR
     7 THE RANGES, THE USER MAY SPECIFY ANOTHER'/' OPTION
     8 WHICH ALLOWS THE PROGRAM TO AUTOMATICALLY FORM'/' THE
     9 BREAKDOWNS BY INDIVIDUAL VALUES OF THE')
632   FORMAT(' BREAKDOWN VARIABLE.  IF MORE VALUES'/
     2 ' EXIST THAN GROUPS, THE PROGRAM AUTOMATICALLY'/' FORMS
     3 THE RANGES FOR BREAKDOWN'/ ' OPTIONS ARE:'/
     5 ' "BREAK" - FORM GROUPINGS BASED ON VALUES OF BREAKDOWN
     6 VARIABLES.'/' "GROUP" - SPECIFY MAXIMUM NUMBER OF GROUPS
     7 (PRESET TO 20) MAXIMUM OF 20.'/' "AUTO" - AUTOMATIC BREAKDOWN
     8 (SPECIFY THIS ONLY WHEN ASKED FOR RANGES)'/' "DISCR" - BOTH',
     8' GROUPS ARE AUTO (SPECIFY NOW)'/' "HEADR" -',
     9' OMIT SIZE, MEANS, AND STD. DEV. REPORT'/' "RANGE" - TYPE OUT',
     1' RANGES FOR BREAKDOWNS IF AUTOMATICALLY FORMED'/
     2'0IF NO OPTIONS ARE DESIRED TYPE A RETURN')
3     FORMAT (5(A5,1X))
      IF(OPT(1).EQ.'!') RETURN
      DO 4 I=1,5
      IF (OPT(I).EQ.'    ') GO TO 10
      IF (OPT(I).NE.'HELP') GO TO 6
      WRITE (IDLG,2)
      WRITE(IDLG,632)
      GO TO 5
6     IF (OPT(I).NE.'BREAK') GO TO 7
      OPTD(1)=1
      GO TO 4
7     IF (OPT(I).NE.'GROUP') GO TO 90
      OPTD(2)=1
      GO TO 4
90    IF(OPT(I).NE.'HEADR') GO TO 308
      OPTD(3)=1
      GO TO 4
308   IF(OPT(I).NE.'DISCR') GO TO 508
      OPTD(4)=1
      GO TO 4
508   IF(OPT(I).NE.'RANGE') GO TO 509
      OPTD(5)=1
      GO TO 4
509   IF(OPT(I).NE.'AUTO') GO TO 512
510   WRITE(IDLG,511)
511   FORMAT(' "AUTO" ONLY SPECIFIED WHEN RANGES ARE CALLED FOR')
      GO TO 5
512   IF(OPT(I).EQ.'AUTO,') GO TO510
8     WRITE (IDLG,9) OPT(I)
9     FORMAT ('0OPTION ',A5,' DOES NOT EXIST TRY AGAIN')
      GO TO 5
4     CONTINUE
10    IF (OPTD(1).EQ.1) GO TO 30
C
C     CELLS INDIVIDUAL VARIABLES.
C
11    IF(ICC.NE.2) WRITE (IDLG,12)
12    FORMAT ('0HOW MANY CELLS IN FACTOR 1? ',$)
      READ (ICC,13) NGRPA
13    FORMAT(I)
      IF (NGRPA.LE.20) GO TO 15
      WRITE (IDLG,14)
14    FORMAT ('0MAXIMUM OF 20')
      GO TO 11
15    IF(ICC.NE.2) WRITE (IDLG,16)
16    FORMAT ('0HOW MANY CELLS IN FACTOR 2? ',$)
      READ (ICC,13) NGRPB
      IF (NGRPB.LE.20) GO TO 18
      WRITE (IDLG,14)
      GO TO 15
18    IF(ICC.NE.2) WRITE (IDLG,19)
19    FORMAT ('0TYPE IN EACH VARIABLE AFTER THE'/
     1' CORRESPONDING  LEVEL, "EMPTY"-INDICATES EMPTY CELL')
      XT=NC
      L=-1
      DO 26 I=1,NGRPA
      DO 26 J=1,NGRPB
22    IF(ICC.NE.2) WRITE (IDLG,21)I,J
21    FORMAT ('0CELL('I2,',',I2,')? ',$)
      IRET=-99
      CALL ALPHA(ICELL,1,K,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IHELP.EQ.1) GO TO 22
      IF(IERR.EQ.1) GO TO 22
      IF(ICELL.GE.0) GO TO 441
      L=L-1
      IF(L.GE.-NV) GO TO 28
445   WRITE(IDLG,27)
27    FORMAT(' TOO MANY CELLS FOR THE NUMBER OF VARIABLES')
      GO TO 11
28    ICELL=L
      GO TO 24
441   K=(I-1)*NGRPB+J-1
      IF(K.LT.1) GO TO 24
      DO 442 II=1,K
      IF(IVA(II).LE.0) GO TO 442
      IF(IVA(II).NE.ICELL) GO TO 442
      WRITE(IDLG,443)
443   FORMAT(' THIS VARIABLE USED PREVIOUSLY')
      GO TO 22
442   CONTINUE
24    K=(I-1)*NGRPB+J
26    IVA(K)=ICELL
      NGPALL=NGRPA*NGRPB
      K=0
      DO 421 I=1,NGPALL
      IF(IVA(I).EQ.0) GOTO 421
      K=K+1
421   CONTINUE
      IF(K.GT.NV) GO TO 445
      GOTO 408
400   J=NGPALL
404   IF(IVA(J).GT.0) GO TO 405
      IVA(J)=IVA(J)-1
      IF(IVA(J).GE.-NV) GO TO 406
405   J=J-1
      IF(J.GE.1) GO TO 404
      RETURN
406   K=IVA(J)
      IF(J.EQ.NGPALL) GO TO 408
      DO 407 I=J+1,NGPALL
      IF(IVA(I).GT.0) GO TO 407
      K=K-1
      IF(K.LT.-NV) GO TO 405
      IVA(I)=K
407   CONTINUE
408   DO 403 I=1,NGPALL-1
      DO 403 J=I+1,NGPALL
      IF((IVA(I).EQ.0).OR.(IVA(J).EQ.0)) GO TO 403
      IF(IVA(I).LT.0) GO TO 409
      IF(IVA(J).GT.0) GOTO 410
      GO TO 411
409   IF(IVA(J).LT.0) GO TO 410
411   IF(IVA(J).EQ.-IVA(I)) GO TO 400
      GO TO 403
410   IF(IVA(J).EQ.IVA(I)) GO TO 400
403   CONTINUE
      DO 88 I=1,20
      AM(I)=0
      Y1(I)=0
      S1(I)=0
      C(I)=0
      B1(I)=0
      S2(I)=0
      Y2(I)=0
      DO 88 J=1,20
      T(I,J)=0
      YT(I,J)=0
88    YS(I,J)=0
      RAB=0
      SST=0
      SS=0
      G=0
      SSC=0
      SSAB=0
      AN=0
      SSBA=0
      TW=0
      TC=0
      DO 20 I=1,NGRPA
      DO 25 J=1,NGRPB
      K=(I-1)*NGRPB+J
      ICELL=IVA(K)
      IF(ICELL.EQ.0) GO TO 20
      IF(ICELL.LT.0) ICELL=-ICELL
23    TC=TC+1.
      TW=TW+XT-1.
      YT(I,J)=VMN(ICELL)*XT
      T(I,J)=XT
      YS(I,J)=((STD(ICELL)**2)*(XT*(XT-1.))+YT(I,J)**2)/XT
      Y1(I)=Y1(I)+YT(I,J)
      S1(I)=S1(I)+T(I,J)
      SSC=SSC+(YT(I,J)*YT(I,J))/T(I,J)
      SST=SST+YS(I,J)
25    CONTINUE
      IF(S1(I).EQ.0) GO TO 20
      SSAB=SSAB+(Y1(I)**2)/S1(I)
      SS=SS+S1(I)
20    G=G+Y1(I)
      GO TO 70
C
C     CELLS ON INDIV. VAR. BY BREAK-DOWN.
C
30    IF (OPTD(2).NE.1) GO TO 36
31    IF(ICC.NE.2) WRITE (IDLG,32)
32    FORMAT ('0MAXIMUM NUMBER OF GROUPS FOR LEVEL 1? ',$)
      READ (ICC,13) NGRPA
      IF ((NGRPA.GT.1).AND.(NGRPA.LE.20)) GO TO 34
      WRITE (IDLG,33)
33    FORMAT ('0NO. OF GROUPS SPECIFIED WAS EITHER LESS THAN 1'/
     1' OR GREATTER THAN THE MAXIMUM (20)')
      GO TO 31
34    IF(ICC.NE.2) WRITE (IDLG,35)
35    FORMAT ('0MAXIMUM NUMBER OF GROUPS FOR LEVEL 2? ',$)
      READ (ICC,13) NGRPB
      IF ((NGRPB.GT.1).AND.(NGRPB.LE.20)) GO TO 36
      WRITE (IDLG,33)
      GO TO 34
36    IF(ICC.NE.2) WRITE (IDLG,38)
38    FORMAT ('0WHICH VARIABLES ARE TO BE ANALYSED? ',$)
      CALL ALPHA(ICELLS,40,NZZ,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IHELP.EQ.1) GO TO 36
      IF(IERR.EQ.1) GO TO 36
      DO 302 I=1,NZZ
      IF(ICELLS(I).NE.-1) GO TO 302
      ALL=1
      NZZ=NV
      GO TO 39
302   CONTINUE
39    DO 40 I=1,2
42    IF(ICC.NE.2) WRITE (IDLG,41) I
41    FORMAT ('0WHICH VARIABLE IS BREAKDOWN VARIABLE',I2,'? ',$)
      CALL ALPHA(ICELL,1,J,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IHELP.EQ.1) GO TO 42
      IF(IERR.EQ.1) GO TO 42
      IF(ICELL.GT.0)GO TO 40
      WRITE(IDLG,450)
450   FORMAT(' *, ?, ALL MAY NOT BE USED HERE')
      GO TO 42
40    IBK(I)=ICELL
      DO 45 I=1,2
      IF(OPTD(4).EQ.1) GO TO 50
44    IF(ICC.NE.2) WRITE (IDLG,46) NAMES(IBK(I))
46    FORMAT ('0LIST THE RANGES FOR BREAKDOWN VARIABLE: ',A5/)
      DO 43 J=1,20
57    READ (ICC,47,END=37,ERR=37) HELP
47    FORMAT (A4,70X)
      IF(HELP.EQ.' ') GO TO 37
      IF(HELP.EQ.'!') RETURN
      IF (HELP.NE.'HELP') GO TO 59
      WRITE (IDLG,48)
48    FORMAT ('IF AUTOMATIC BREAK-DOWN BY GROUPS IS DESIRED TYPE "AUTO",
     1'/' OTHERWISE ENTER THE RANGE FOR EACH BREAK-DOWN SMALLER FIRST,
     2  1 PER LINE.'/'EXAMPLE:'/'17,35.2'/'CONTROL Z (^Z) WILL TERMINATE
     3 ENTRY')
      GO TO 44
59    IF (HELP.NE.'AUTO') GO TO 55
300   IF (J.EQ.1) GO TO 50
      WRITE (IDLG,49)
49    FORMAT ('0AUTO CAN ONLY BE SPECIFIED FIRST TIME'/' CONFUSION
     1 START OVER WITH SAME VARIABLE')
      GO TO 44
55    REREAD 56, BRNG(I,J,1),BRNG(I,J,2)
56    FORMAT (2F)
      IF (BRNG(I,J,1).LE.BRNG(I,J,2)) GO TO 43
      WRITE (IDLG,58)
58    FORMAT (' SMALLER FIRST, PLEASE RETYPE'/)
      GO TO 57
43    CONTINUE
      J=J+1
37    NGRPS(I)=J-1
      IF(NGRPS(I).LT.1) RETURN
      GO TO 45
50    XMAX=-10.**11
      XMIN=10.**11
      N=0
      L=IBK(I)
      DO 51 K=1,NC
      IF (XMAX.LT.DATA(K,L))XMAX=DATA(K,L)
      IF (XMIN.GT.DATA(K,L))XMIN=DATA(K,L)
      IF (N.GT.20) GO TO 51
      IF (N.LT.1) GO TO 53
      DO 52 M=1,N
      IF (DATA(K,L).EQ.BRNG(I,M,1)) GO TO 51
52    CONTINUE
53    N=N+1
      IF(N.GT.20) GO TO 51
      BRNG (I,N,1)=DATA(K,L)
      BRNG (I,N,2)=DATA(K,L)
51    CONTINUE
      IF (N.LT.NGRPS(I)) NGRPS(I)=N
      IF (N.EQ.NGRPS(I)) GO TO 89
      XT=NGRPS(I)
      RG=(XMAX-XMIN)/XT
      BRNG(I,1,1)=XMIN
      BRNG(I,1,2)=XMIN+RG
      DO 54 K=2,NGRPS(I)
      BRNG(I,K,1)=BRNG(I,K-1,2)
      BRNG(I,K,2)=BRNG(I,K,1)+RG
54    CONTINUE
      BRNG(I,NGRPS(I),2)=XMAX
89    IF(OPTD(5).NE.1) GO TO 420
      IF(OPTD(4).EQ.1) WRITE(IDLG,301) NAMES(IBK(I))
301   FORMAT('0BREAKDOWN RANGES FOR VARIABLE: ',A5/)
C     *******************************************************
C     SORT BEFORE PUTTING OUT RANGES (BUBBLE SORT BECAUSE OF SIZE)
C
420   J=0
310   J=J+1
      IF(J.EQ.NGRPS(I)) GO TO 312
      TAB=BRNG(I,J+1,1)
      TBA=BRNG(I,J+1,2)
      IF(BRNG(I,J,1).LT.TAB) GO TO 310
      K=J
311   BRNG(I,K+1,1)=BRNG(I,K,1)
      BRNG(I,K+1,2)=BRNG(I,K,2)
      K=K-1
      IF(K.LT.1) GO TO 313
      IF(TAB.LT.BRNG(I,K,1)) GO TO 311
313   BRNG(I,K+1,1)=TAB
      BRNG(I,K+1,2)=TBA
      GO TO 310
C     ***********************************************************
312   IF(OPTD(5).NE.1) GO TO45
      DO 71 J=1,NGRPS(I)
71    WRITE(IDLG,72)(BRNG(I,J,K),K=1,2)
72    FORMAT('+',G10.4,',',G10.4/)
      GO TO 45
45    CONTINUE
C
C     AT THIS POINT RANGES HAVE BEEN DET.
C     BREAK DOWN VARIABLES
C
      DO 2000 KCELL=1,NZZ
      IF(ALL.EQ.1) GO TO 306
      ICELL=ICELLS(KCELL)
      GO TO 305
306   IF((KCELL.EQ.IBK1).OR.(KCELL.EQ.IBK2)) GO TO 2000
      ICELL=KCELL
305   DO 388 I=1,20
      AM(I)=0
      Y1(I)=0
      S1(I)=0
      C(I)=0
      B1(I)=0
      S2(I)=0
      Y2(I)=0
      DO 388 J=1,20
      T(I,J)=0
      YT(I,J)=0
388   YS(I,J)=0
      RAB=0
      SST=0
      SS=0
      G=0
      SSC=0
      SSAB=0
      AN=0
      SSBA=0
      TW=0
      TC=0
      DO 60 I=1,NC
C     ----------------------------------------------------
C     WHICH RANGE IN BREAKDOWN VARIABLE 1 DOES IT BELONG TO
      DO 61 J=1,NGRPA
      IF (DATA(I,IBK1).LT.BRNG(1,J,1)) GO TO 61
      IF (DATA(I,IBK1).GT.BRNG(1,J,2)) GO TO 61
C     -----------------------------------------------------------
C     WHICH RANGE IN BREAKDOWN VARIABLE 2 DOES IT BELONG
62    DO 63 K=1,NGRPB
      IF (DATA(I,IBK2).LT.BRNG(2,K,1)) GO TO 63
      IF (DATA(I,IBK2).GT.BRNG(2,K,2)) GO TO 63
C     --------------------------------------------------------
C     ADD IT INTO THE CORRECT TOTALS
64    YT(J,K)=YT(J,K)+DATA(I,ICELL)
      YS(J,K)=YS(J,K)+DATA(I,ICELL)**2
      T(J,K)=T(J,K)+1.
      GO TO 60
63    CONTINUE
      GO TO 60
61    CONTINUE
60    CONTINUE
92    DO 65 I=1,NGRPA
      DO 66 J=1,NGRPB
      IF(T(I,J).EQ.0) GO TO 66
      TC=TC+1
      TW=TW+T(I,J)-1.
      Y1(I)=Y1(I)+YT(I,J)
      S1(I)=S1(I)+T(I,J)
      SSC=SSC+YT(I,J)**2/T(I,J)
      SST=SST+YS(I,J)
66    CONTINUE
      IF(S1(I).EQ.0) GO TO 65
      SSAB=SSAB+(Y1(I)*Y1(I))/S1(I)
      SS=SS+S1(I)
65    G=G+Y1(I)
      GO TO 70
C
C     BEGIN ANALYSIS OF VARIANCE.
C
70    N1=NGRPA
      N2=NGRPB
      DO 93 J=1,N2
      DO 93 I=1,N1
93    S2(J)=S2(J)+T(I,J)
      SWTCH=0
      DO 73 I=1,N1
      IF(S1(I).NE.0) GO TO 73
      IF(SWTCH.EQ.1) GO TO 75
      WRITE(IDLG,74)
      SWTCH=1
74    FORMAT('OTHE FOLLOWING LEVELS ARE EMPTY - WILL BE ELIMIN.')
75    WRITE(IDLG,76) I
76    FORMAT(' FACTOR 1 LEVEL',I3)
73    CONTINUE
      DO 77 J=1,N2
      IF(S2(J).NE.0) GO TO 77
      IF(SWTCH.EQ.1) GO TO 78
      WRITE(IDLG,74)
      SWTCH=1
78    WRITE(IDLG,79) J
79    FORMAT(' FACTOR 2 LEVEL',I3)
77    CONTINUE
      IF(SWTCH.EQ.0) GO TO 80
81    DO 82 I=1,NGRPA
      IF(S1(I).NE.0) GO TO 82
      NGRPA=NGRPA-1
      IF(I.GT.NGRPA) GO TO 84
      DO 83 K=I,NGRPA
      DO 83 J=1,NGRPB
      YT(K,J)=YT(K+1,J)
      YS(K,J)=YS(K+1,J)
      T(K,J)=T(K+1,J)
      S1(K)=S1(K+1)
83    CONTINUE
      GO TO 81
82    CONTINUE
84    DO 85 J=1,NGRPB
      IF(S2(J).NE.0) GO TO 85
      NGRPB=NGRPB-1
      IF(J.GT.NGRPB) GO TO 87
      DO 86 K=J,NGRPB
      DO 86 I=1,NGRPA
      YT(I,K)=YT(I,K+1)
      YS(I,K)=YS(I,K+1)
      T(I,K)=T(I,K+1)
      S2(K)=S2(K+1)
86    CONTINUE
      GO TO 84
85    CONTINUE
87    DO 91 I=1,20
      Y1(I)=0
      AM(I)=0
      S1(I)=0
      C(I)=0
      B1(I)=0
      S2(I)=0
91    Y2(I)=0
      RAB=0
      SST=0
      SS=0
      G=0
      SSC=0
      SSAB=0
      AN=0
      SSBA=0
      TW=0
      TC=0
      GO TO 92
80    DO 100 J=1,N2
      DO 101 I=1,N1
      IF(T(I,J).EQ.0) GO TO 101
102   Y2(J)=Y2(J)+YT(I,J)
      YS(I,J)=YS(I,J)-(YT(I,J)*YT(I,J))/T(I,J)
      IF(T(I,J).EQ.1.) GO TO 101
106   YS(I,J)=YS(I,J)/(T(I,J)-1.)
      IF(YS(I,J).LT.0) YS(I,J)=0.
      YS(I,J)=SQRT(YS(I,J))
      YT(I,J)=YT(I,J)/T(I,J)
101   CONTINUE
      SSBA=SSBA+(Y2(J)*Y2(J))/S2(J)
100   CONTINUE
      G1=(G*G)/SS
      SST=SST-G1
      SSC=SSC-G1
      SSAB=SSAB-G1
      SSBA=SSBA-G1
      SSW=SST-SSC
      NP=0
      NP=NP+1
      IF(TC.EQ.1.) GO TO 149
104   SMC=SSC/(TC-1.)
      NP=NP+1
      ITYPE=0
      IF(TW.EQ.0) ITYPE=1
105   IF(ITYPE.EQ.1) GO TO 505
      SMW=SSW/TW
      NP=NP+1
      IF(SMW.EQ.0) GO TO 149
103   F=SMC/SMW
505   IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,321)
321   FORMAT('0',15X,'*****2-WAY ANOVA*****')
      LINES=4
      IF(OPTD(1).NE.1) GO TO 323
      WRITE(IOUT,320)NAMES(ICELL),NAMES(IBK1),NAMES(IBK2)
320   FORMAT(' ANALYSIS RUN ON VARIABLE: ',A5,' WITH CELLS DETERMINED'/
     1' BY BREAKDOWNS ON VARIABLE: ',A5,' AND VARIABLE: ',A5)
      LINES=LINES+2
323   IF(OPTD(3).EQ.1) GO TO 155
      LINES=LINES+3+(N2+ISQ-1)/ISQ
324   WRITE (IOUT,107)
107   FORMAT('0FACTOR'/2X,'ONE',25X,'FACTOR TWO')
      DO 322 K=1,N2,ISQ
      NEND=K+ISQ-1
      IF(NEND.GT.N2) NEND=N2
322   WRITE(IOUT,108)(I,I=K,NEND)
108   FORMAT(' LEVEL',5X,15(4X,I3,1X))
      DO 116 I=1,N1
      IF(IOUT.NE.21) GO TO 325
      LINES=LINES+3*((N1+ISQ-1)/ISQ)+1
      IF(LINES.LE.LINPP) GO TO 325
      CALL PRNTHD
      WRITE(IOUT,107)
      DO 326 K=1,N2,ISQ
      NEND=K+ISQ-1
      IF(NEND.GT.N2) NEND=N2
326   WRITE(IOUT,108)(IK,IK=K,NEND)
      LINES=5+((N2+ISQ-1)/ISQ)*4+1
325   DO 116 K=1,N2,ISQ
      NEND=K+ISQ-1
      IF(NEND.GT.N2) NEND=N2
      IF(K.EQ.1) GO TO 112
      WRITE(IOUT,110)(T(I,J),J=K,NEND)
110   FORMAT('0',5X,'N',4X,15F8.2)
      GO TO 111
112   WRITE (IOUT,113)I,(T(I,J),J=K,NEND)
113   FORMAT('0',I3,2X,'N',4X,15F8.2)
111   WRITE (IOUT,114)(YT(I,J),J=K,NEND)
114   FORMAT(6X,'MEAN',1X,15F8.2)
115   WRITE (IOUT,154)(YS(I,J),J=K,NEND)
154   FORMAT(6X,'STDV',1X,15F8.2)
116   CONTINUE
      IF(ITYPE.EQ.1) GO TO 335
155   IF(IOUT.NE.21) GO TO 327
      LINES=LINES+9
      IF(LINES.LE.LINPP) GO TO 327
      CALL PRNTHD
      LINES=11
327   WRITE (IOUT,117)
117   FORMAT(//24X,'PRELIMINARY ANOVA')
      WRITE (IOUT,118)
118   FORMAT(6X,'SOURCE',11X,'DF',8X,'SS',8X,'MS',8X,'F',4X,'PROB')
      KK=TC-1.
      NDG=TW
      PROB=FISHER(KK,NDG,F)
      WRITE (IOUT,119)KK,SSC,SMC,F,PROB
119   FORMAT(6X'CELLS',9X,I5,3F10.2,2X,F6.4)
      KK=N1-1
      I=1
      J=2
      WRITE (IOUT,120)I,J,KK,SSAB
120   FORMAT(4X,I1,' IGNORING',I2,4X,I5,F10.2)
      KK=N2-1
      WRITE (IOUT,120)J,I,KK,SSBA
      KK=TW
      WRITE (IOUT,121)KK,SSW,SMW
121   FORMAT(6X,'WITHIN',8X,I5,2F10.2)
      KK=SS-1.
      WRITE (IOUT,122)KK,SST
122   FORMAT(6X,'TOTAL',9X,I5,F10.2)
335   DO 125 I=1,N2
      DO 125 J=1,N2
      YT(I,J)=0.
      IF(I-J)124,123,124
123   YS(I,J)=1.
      GO TO 125
124   YS(I,J)=0.
125    CONTINUE
      DO 128 I=1,N2
      DO 128 J=1,N2
      DO 128 K=1,N1
      IF(S1(K))128,126,128
126   TYPE 127,K,S1(K)
127   FORMAT(3HS1(,I3,2H)=,F6.0)
128   YT(I,J)=YT(I,J)+(T(K,I)*T(K,J))/S1(K)
      DO 129 I=1,N2
129   YT(I,I)=YT(I,I)-S2(I)
      DO 133 J=1,N2
130   DO 132 K=1,N1
      IF(S1(K))132,131,132
131   TYPE 127,K,S1(K)
132   C(J)=C(J)+(T(K,J)*Y1(K))/S1(K)
      YT(N2,J)=1.
133   C(J)=C(J)-Y2(J)
      C(N2)=0.
C
C     *************************************************************
C     INVERSE BY LINEAR ROW TRANSFORMATION
C
      DO 202 I=1,N2
      DO 201 J=1,N2
201   YS(I,J)=0
202   YS(I,I)=1.
      DO 210 I=1,N2
      IF(((YT(I,I)+100.)-100.).NE.0.0) GO TO 220
      IF(I.EQ.N2) GO TO 260
      DO 211 J=I+1,N2
      IF(((YT(J,I)+100.)-100.).NE.0.0) GO TO 212
211   CONTINUE
260   WRITE(IDLG,261)
261   FORMAT('0INVERSE DOES NOT EXIST')
      GO TO 1000
212   DO 213 K=1,N2
      YT(I,K)=YT(I,K)+YT(J,K)
213   YS(I,K)=YS(I,K)+YS(J,K)
220   GOB=YT(I,I)
      DO 221 J=1,N2
      YT(I,J)=YT(I,J)/GOB
221   YS(I,J)=YS(I,J)/GOB
      DO 230 L=1,N2
      IF(L.EQ.I) GO TO 230
      GOB=YT(L,I)
      DO 231 J=1,N2
      YT(L,J)=YT(L,J)-GOB*YT(I,J)
231   YS(L,J)=YS(L,J)-GOB*YS(I,J)
230   CONTINUE
210   CONTINUE
C     ***********************************************************
C
270   DO 134 J=1,N2
      DO 134 K=1,N2
134   B1(J)=B1(J)+YS(J,K)*C(K)
      DO 138 I=1,N1
      DO 135 J=1,N2
135   AM(I)=AM(I)+T(I,J)*B1(J)
      AM(I)=Y1(I)-AM(I)
      IF(S1(I))137,136,137
136   TYPE 127,I,S1(I)
137   AM(I)=AM(I)/S1(I)
138   AN=AN+AM(I)
      P1=N1
      AN=AN/P1
      DO 139 I=1,N1
      A=AM(I)-AN
139   RAB=RAB+A*Y1(I)
      DO 140 J=1,N2
140   RAB=RAB+B1(J)*Y2(J)
      RAB=RAB+AN*G-G1
      SI=SSC-RAB
      P2=N2
      TI=TC-P1-P2+1.
      MP=0
      MP=MP+1
      IF(TI)141,151,141
141   SMI=SI/TI
      TAB=RAB-SSBA
      MP=MP+1
      IF(P1-1.)142,151,142
142   TMAB=TAB/(P1-1.)
      TBA=RAB-SSAB
      MP=MP+1
      IF(P2-1.)143,151,143
143   TMBA=TBA/(P2-1.)
      MP=MP+1
      NDG=TW
      IF(ITYPE.EQ.1) NDG=TI
      IF (ITYPE.EQ.1) SMW=SMI
      IF(SMW)144,151,144
144   FI=SMI/SMW
      FA=TMAB/SMW
      FB=TMBA/SMW
      IF(IOUT.NE.21) GO TO 328
      LINES=LINES+10
      IF(LINES.LE.LINPP) GO TO 328
      CALL PRNTHD
      LINES=10
328   WRITE (IOUT,145)
145   FORMAT(//29X11HFINAL ANOVA)
      WRITE (IOUT,118)
      KK=TC-1.
      WRITE (IOUT,119)KK,SSC
      KK=P1-1.
      I=1
      J=2
      PROB=FISHER(KK,NDG,FA)
      WRITE (IOUT,146)I,J,KK,TAB,TMAB,FA,PROB
146   FORMAT(3XI1,12H ELIMINATING,I2,2X,I5,3F10.2,2X,F6.4)
      KK=P2-1.
      PROB=FISHER(KK,NDG,FB)
      WRITE (IOUT,147)J,I,KK,TBA,TMBA,FB,PROB
147   FORMAT(3XI1,12H ELIMINATING,I2,2X,I5,3F10.2,2X,F6.4)
      KK=TI
      IF(ITYPE.EQ.0) PROB=FISHER(KK,NDG,FI)
      IF(ITYPE.EQ.0)WRITE (IOUT,148)I,J,KK,SI,SMI,FI,PROB
      IF(ITYPE.EQ.1) WRITE(IOUT,148)I,J,KK,SI,SMI
148   FORMAT(6XI1,3H BY,I2,8X,I5,3F10.2,2X,F6.4)
      KK=TW
      IF(ITYPE.EQ.0)WRITE (IOUT,121)KK,SSW,SMW
      KK=SS-1.
      WRITE (IOUT,122)KK
      GO TO 1000
C     VALUES OF NP MAY BE 1,2,3
C     1 IMPLIES TC-1=0,2 IMPLIES TW=0,3 IMPLIES SMW=0
149   TYPE 150,NP
150   FORMAT(3HNP=,I2)
      GO TO(104,105,103),NP
C     VALUES OF MP MAY BE 1,2,3,4
C     1 IMPLIES TI=0,2 IMPLIES P1-1=0,3 IMPLIES P2-1=0,4 IMPLIES SMW=0
151   TYPE152,MP
152   FORMAT(3HMP=,I2)
      GO TO(141,142,143,144),MP
1000  IF(OPTD(1).NE.1) GO TO 400
2000  CONTINUE
      RETURN
      END
C                                                     *** STAT PACK ***
C     SUBROUTINE FOR FRIEDMAN TWO-WAY ANALYSES OF VARIANCE
C     CALLING SEQUENCE: CALL FRIED(NV,NC,MV,MC,DATA,SUM,RANK,NAMES)
C     WHERE: NV - NUMBER OF VARIABLES USED
C            NC - NUMBER OF OBSERVATIONS USED
C            MV - MAXIMUM NUMBER OF VARIABLES ALLOWED
C            MC - MAXIMUM NUMBER OF OBSERVATIONS ALLOWED
C            DATA - MATRIX CONTAINING DATA
C            SUM - EXTRA VECTOR AT LEAST NV LONG
C            RANK - EXTRA VECTOR AT LEAST NV LONG
C            NAMES - VECTOR CONTAINING VARIABLE NAMES
C
C     ROUTINE SUGGESTED TO BE INCLUDED IN STAT PACK BY LONNIE
C     HANNAFORD (SPECIAL EDUCATION) AND ULDIS SMIDCHENS (EDUCATION).
C     SOURCE IS NON PARAMETRIC STATISTICS BY SIEGEL PAGES 166-172.
C
      SUBROUTINE FRIED(NV,NC,MV,MC,DATA,SUM,RANK,NAMES)
      DIMENSION DATA(MC,MV),SUM(1),RANK(1),IV(20),IVB(20),NAMES(1)
      COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON /EXTRA/ HEDR(70),NSZ
1     IF(ICC.NE.2) WRITE(IDLG,2)
2     FORMAT(' WHICH VARIABLES? '$)
      CALL ALPHA(IV,20,N,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IERR.EQ.1) GO TO 1
      IF(IHELP.NE.1) GO TO 4
      WRITE(IDLG,3)
3     FORMAT('0ENTER UP TO 20 VARIABLES FOR THE FRIEDMAN TWO-WAY'/
     1' ANALYSIS OF VARIANCE.  VARIABLE NAMES OR NUMBERS'/
     2' MAY BE USED TO INDICATED THE VARIABLES.  RANGES OF'/
     3' VARIABLES MAY BE ENTERED BY TYPING THE EXTREMES OF THE RANGE'/
     4' SEPARATED BY A MINUS.  THE ASTERIK MAY BE USED IN ONE'/
     5' OR MORE POSITIONS TO INDICATED ALL POSSIBLE COMBINATIONS.')
      GO TO 1
4     IF(N.LT.2) RETURN
      DO 5 I=1,N-1
      IF(IV(I).LE.0) GO TO 5
      DO 6 J=I+1,N
      IF(IV(I).NE.IV(J)) GO TO 6
      WRITE(IDLG,7) NAMES(IV(J))
7     FORMAT(' SAME VARIABLE (',A5,') LISTED TWICE')
      GO TO 1
6     CONTINUE
5     CONTINUE
C
C     OK NOW INSERT FOR * AND ALL.
C
      K=1
      DO 8 I=1,N
      IVB(I)=IV(I)
      IF(IV(I).GT.0) GO TO 8
      IVB(I)=K
      K=K+1
8     CONTINUE
      GO TO 18
C
C     RETURN HERE TO PICK UP NEXT SET OF VARIABLES
C
10    J=N
14    IF(IV(J).GT.0) GO TO 15
      IVB(J)=IVB(J)+1
      IF(IVB(J).LE.NV) GO TO 16
15    J=J-1
      IF(J.GE.1) GO TO 14
      RETURN
16    K=IVB(J)
      IF(J.EQ.N) GO TO 18
      DO 17 I=J+1,N
      IF(IV(I).GT.0) GO TO 17
      K=K+1
      IF(K.GT.NV) GO TO 15
      IVB(I)=K
17    CONTINUE
18    DO 13 I=1,N-1
      DO 13 K=I+1,N
      IF(IVB(I).EQ.IVB(K)) GO TO 14
13    CONTINUE
C
C     BEGIN THE FRIEDMAN ANALYSIS OF VARIANCE
C
      DO 19 I=1,N
19    SUM(I)=0
      DO 20 J=1,NC
      DO 21 I=1,N
21    RANK(I)=DATA(J,IVB(I))
C
C     RANKING ( ONLY 20 MAX, SO SORT RANK METHOD NOT USED )
C
      DO 22 I=1,N
      SAME=1.
      ALESS=0
      DO 23 K=1,N
      IF(RANK(K).GT.RANK(I)) GO TO 23
      IF(RANK(K).LT.RANK(I)) GO TO 24
      SAME=SAME+1
      GO TO 23
24    ALESS=ALESS+1.
23    CONTINUE
22    SUM(I)=SUM(I)+ALESS+SAME/2.
20    CONTINUE
C
C     RANKING COMPLETE NOW CALCULATE THE CHI SQUARE
C
      IF(IOUT.NE.21) WRITE(IOUT,5566) (HEDR(J),J=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,31)
31    FORMAT('0',8X,'**** FRIEDMAN TWO-WAY ANALYSIS OF VARIANCE ****'/
     1'0VARIABLE  SUM OF RANKS'/' ========  ============')
      LINES=7
      TOTSQ=0
      DO 30 I=1,N
      IF(IOUT.NE.21) GO TO 36
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 36
      CALL PRNTHD
      WRITE(IOUT,31)
      LINES=8
36    WRITE(IOUT,32) NAMES(IVB(I)),SUM(I)
32    FORMAT(2X,A5,5X,F10.2)
30    TOTSQ=TOTSQ+SUM(I)**2
      FRIED=(12./(NC*N*(N+1.)))*TOTSQ-3.*NC*(N+1.)
      IDF=N-1
      FD=FRIED/IDF
      PROB=FISHER(IDF,1000,FD)
      IF(IOUT.NE.21) GO TO 37
      LINES=LINES+3
      IF(LINES.LE.LINPP) GO TO 37
      CALL PRNTHD
37    WRITE(IOUT,35) FRIED,IDF,PROB
35    FORMAT('0CHI SQUARE=',F12.3,', WITH ',I2,' DEGREES OF FREEDOM'/
     1' HAVING A PROBABILITY OF ',F5.2)
      GO TO 10
      END
C                                                     *** STAT PACK ***
C     SUBROUTINE FOR SIGN TEST.
C     CALLING SEQUENCE: CALL SIGNT(NV,NC,MV,MC,DATA,NAMES)
C     WHERE: NV - NUMBER OF VARIABLES USED.
C            NC - NUMBER OF OBSERVATIONS USED.
C            MV - MAXIMUM NUMBER OF VARIABLES ALLOWED
C            MC - MAXIMUM NUMBER OF OBSERVATIONS ALLOWED
C            DATA - MATRIX FOR DATA
C            NAMES - VECTOR CONTAINING VARIABLE NAMES.
C
C     ROUTINE SUGGESTED BY LONNIE HANNAFORD (SPECIAL EDUCATION) AND
C     ULDIS SMIDCHENS (TEACHER EDUCATION).  SOURCE IS NON-
C     PARAMETRIC STATISTICS BY SIEGEL PAGES 68-75.  THE
C     BINOMIAL EXPANSION COULD BE USED FOR N>25, HOWEVER IT IS
C     SLOW FOR LARGE N.  WORKS BY ATTEMPTING TO NORMALIZE NUMBER
C     AS IT IS BEING CALCULATED ABOUT THE VALUE 1.
C
      SUBROUTINE SIGNT(NV,NC,MV,MC,DATA,NAMES)
      DIMENSION DATA(MC,MV),NAMES(1)
      COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON /EXTRA/ HEDR(70),NSZ
C     EQUIVALENCE FOR BINARY EXPANSION USED TO FIND POWER OF 2
      EQUIVALENCE (ISUM,SUM)
      IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,1)
1     FORMAT('0',20X,'**** SIGN TEST ****'/'0',35X,'FREQUENCY OF'/
     119X,'NUMBER OF',9X,'SIGN WHICH'/19X,'PAIRS WITH',6X,
     2'OCCURED LEAST',7X,'Z',10X,'PROB'/' VAR',4X,'VAR',4X,
     3'SIGNED DIFFERENCE',7X,'OFTEN',8X,'(N>25)',6X,'(N<=25)')
      LINES=9
      DO 2 I=1,NV-1
      DO 2 J=I+1,NV
      N=0
      NMIN=0
      DO 3 K=1,NC
      IF(DATA(K,I)-DATA(K,J)) 4,3,5
4     NMIN=NMIN+1
5     N=N+1
3     CONTINUE
      IF(NMIN.GT.(N-NMIN)) NMIN=N-NMIN
      IF(N.GT.25) GO TO 20
C
C     BINOMIAL EXPANSION (N<=25)
C
      XSUM=0
      IF(N.LT.127) XSUM=1./2.**N
      IF(NMIN.LT.1) GO TO 13
      DO 11 L=1,NMIN
      SUM=1.
      IN=N-L+1
      XONE=1.
      NTWOS=N
      DO 12 K=IN,N
      SUM=SUM*K/XONE
      IF(NTWOS.LT.1) GO TO 12
      IF(SUM.LT.2) GO TO 12
      JJ=(ISUM.AND."377777777777)/2**27-129
      NDIV=JJ
      IF(JJ.GT.NTWOS) NDIV=NTWOS
      SUM=SUM/2.**NDIV
      NTWOS=NTWOS-NDIV
12    XONE=XONE+1.
      IF(NTWOS.GT.127) GO TO 11
      IF(NTWOS.GT.0) SUM=SUM/2.**NTWOS
      XSUM=XSUM+SUM
11    CONTINUE
13    IF(IOUT.NE.21) GO TO 15
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 15
      CALL PRNTHD
      WRITE(IOUT,1)
      LINES=10
15    WRITE(IOUT,14) NAMES(I),NAMES(J),N,NMIN,XSUM
14    FORMAT(1X,A5,2X,A5,8X,I5,13X,I5,20X,F7.4)
      GO TO 2
C
C     BINOMIAL DISTRIBUTION GIVE Z VALUE
C
20    A=N
      XMN=A/2.
      XSTD=SQRT(A)/2.
      Z=N-XMN
      IF(N.LT.XMN) Z=Z+.5
      Z=Z/XSTD
      IF(IOUT.NE.21) GO TO 22
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 22
      CALL PRNTHD
      WRITE(IOUT,1)
      LINES=10
22    WRITE(IOUT,21) NAMES(I),NAMES(J),N,NMIN,Z
21    FORMAT(1X,A5,2X,A5,8X,I5,13X,I5,7X,G12.3)
2     CONTINUE
      RETURN
      END

