C	WESTERN MICHIGAN UNIVERSITY
C	TRIAOV.F4 (FILENAME ON LIBRARY DECTAPE)
C	TRIAOV, 1.9.3 (CALLING NAME, SUBLST NO.)
C	THREE-WAY ANALYSIS OF VARIANCE
C	PROGRAMMED BY SAM ANEMA AT WMU, LATER MODIFIED BY R.R. BARR,
C	 M.T. O'BRYAN, CONSULTANT - DR. M. STOLINE
C	LIBRARY DECTAPE PROGS. USED:  USAGE.MAC
C	APLIB PROGS. USED:  IO, GETFOR, CHIPRB, FISHER, NORMCV
C	FORWMU PROGS. USED:  TTYPTY, MINSQV, ALLCOR, PRINTS, DEVCHG,
C	 EXISTS
C	INTERNAL SUBR. USED:  MNL, BARTLT, FIT, PRINT, WEIGH, UNWEI,
C	 PROP
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C---------------S(1) USED IN CALL MNL, LX(3) CONTAINS NO. OF LEVELS
C--------------- ON 3 FACTORS.  LNAME(2) IS APPARENTLY NOT USED
      DIMENSION LX(3),S(1),LNAME(2)
C
C     READ IN PARAMETERS
C
99    TYPE 100
100   FORMAT('1 ---WMU 3-WAY ANALYSIS OF VARIANCE --'///)
C	CALL USAGE('TRIAOV')
C---------------TTYPTY RETURNS ZERO - TTY, MINUS ONE - BATCH JOB
	CALL TTYPTY(ICODE)
C---------------ODV, IDV ARE RETURNED.  OTHER ARGS. ARE INPUT.
C--------------- 1 MEANS OUTPUT? PRINTS.  0 MEANS INPUT? PRINTS.
	CALL IO(6,ODV,-1,-4,1,ICODE)
801	CALL IO(4,IDV,-1,-4,0,ICODE)
98    TYPE 101
101   FORMAT (' ENTER FACTOR NAMES AND NO. OF LEVELS'/)
90    ACCEPT 102,INAME,IMAX,JNAME,JMAX
102   FORMAT (A1,1X,I)
      ACCEPT 108,KNAME,KMAX,NTIMES
108   FORMAT(A1,1X,2I)
      IF(IMAX.LT.2.OR.JMAX.LT.2.OR.KMAX.LT.2)GO TO 60
      GO TO 63
60    TYPE 64
64    FORMAT(' THE PROGRAM REQUIRES AT LEAST 2 LEVELS ON EACH 
     1FACTOR.'/' REENTER FACTOR NAMES AND NO. OF LEVELS.'/)
	IF(ICODE.NE.0)CALL EXIT
      GO TO 90
61    TYPE 62
62    FORMAT(' THE PROBLEM IS TOO LARGE.'/' REENTER FACTOR NAMES AND 
     1LEVELS.'/)
	IF(ICODE.NE.0)CALL EXIT
      GO TO 90
63    LX(1)=IMAX
      LX(2)=JMAX
      LX(3)=KMAX
      DO 10 I=2,3
      DO 11 J=1,I-1
      IF(LX(I).GE.LX(J))GO TO 11
      ITM=LX(I)
      LX(I)=LX(J)
      LX(J)=ITM
11    CONTINUE
10    CONTINUE
      MX=IMAX*JMAX*KMAX
      MAX=2*(MX+IMAX*JMAX+IMAX*KMAX+JMAX*KMAX)+10*LX(3)+(LX(3)+1)*
     1	(LX(3)+1)+2*LX(3)*LX(2)+2+IMAX+JMAX+KMAX+MX
      CALL ALLCOR(MAX,IERR,I1,S)
      IF(IERR.NE.0)GO TO 61
      I2=I1+MX
      I3=I2+MX
      I4=I3+IMAX*JMAX
      I5=I4+IMAX*KMAX
      I6=I5+LX(3)+1
      I7=I6+LX(3)+1
      I8=I7+JMAX*KMAX
      I9=I8+IMAX
      I10=I9+JMAX
      I11=I10+KMAX
      I12=I11+IMAX*JMAX
      I13=I12+IMAX*KMAX
      I14=I13+JMAX*KMAX
      I15=I14+LX(3)
      I16=I15+LX(3)
      I17=I16+LX(3)
      I18=I17+(LX(3)+1)*(LX(3)+1)
      I19=I18+LX(3)
      I20=I19+LX(3)
      I21=I20+LX(3)*LX(2)
      I22=I21+LX(3)*LX(2)
      I23=I22+LX(3)
      I24=I23+LX(3)
	I25=I24+LX(3)
      CALL MNL(INAME,JNAME,KNAME,
     1        IMAX,JMAX,KMAX,NTIMES,S(I1),S(I2),S(I3),S(I4),S(I5),
     1S(I6),S(I7),S(I8),S(I9),S(I10),S(I11),S(I12),S(I13),S(I14),S(I15)
     2,S(I16),S(I17),S(I18),S(I19),S(I20),S(I21),S(I22),S(I23),S(I24),
     3	S(I25),IDV)
	TYPE 802
802	FORMAT(//)
	GO TO 801
      END
C---------------INAME, JNAME, KNAME, IMAX, JMAX, KMAX, NTIMES
C--------------- ARE INPUT.  OTHER ARGS. ARE SPACES 
C--------------- RESERVED BY DYN. ALLOC.
      SUBROUTINE MNL(INAME,JNAME,KNAME,
     1               IMAX,JMAX,KMAX,NTIMES,NIJK,XIJK,XIJ,XIK,MC,MR,XJK,
     1XI,XJ,XK,NIJ,NIK,NJK,NI,NJ,NK,D,YB,XD,DF,FF,C,B,A,XXIJK,IDV)
      DIMENSION NIJK(1),XIJK(1),XIJ(1),XIK(1),MC(1),MR(1),XJK(1),XI(1),
     1XJ(1),XK(1),ID(5,50),FMT(48),X(50),NIJ(1),NIK(1),NJK(1),NI(1),
     2NJ(1),NK(1),D(1),YB(1),XD(1),DF(1),FF(1),C(1),B(1),A(1),IAST(19)
     3,IENT(17),XXIJK(1)
      DATA IAST(18),IAST(19)/',0,0,','0,0,0'/
63    IF(NTIMES.EQ.0)NTIMES=1
      TYPE 103
103   FORMAT (' ENTER ID IF DESIRED ELSE RETURN'/)
      DO 37 J=1,NTIMES
37    ACCEPT 104, (ID(I,J),I=1,5)
104   FORMAT (16A5)
      TYPE 1010
1010  FORMAT(' MEANS? ',$)
      ACCEPT 1011,IFM
1011  FORMAT(A3)
C---------------FMT RETURNS USER FORMAT.  ISTD=1 IF STANDARD
C--------------- FORMAT IS REQUESTED, =0 IF FORMAT TO BE USED IS IN FMT.
      CALL GETFOR(-1,-4,FMT,ISTD,48,4)
96    TYPE 2191
2191   FORMAT(' WHICH METHOD OF DATA ENTRY? (1 OR 2) ',$)
      ACCEPT 2192,NWHICH
2192  FORMAT(I)
	IF(IDV.NE.'TTY')GO TO 97
      TYPE 95
95    FORMAT(' ENTER DATA'/)
	GO TO 94
97	TYPE 93
93	FORMAT(' DATA IS BEING READ'/)
C
C     INITIALIZE
C
94      NTI=0
   80 NTI=NTI+1
      N=0
      MX= IMAX*JMAX*KMAX
      DO 200 I=1,MX
      NIJK(I)=0
	XXIJK(I)=0
200   XIJK(I)=0.0
      MX= IMAX*JMAX
      DO 201 I=1,MX
      NIJ(I)=0
201   XIJ(I)=0.0
      MX=JMAX*KMAX
      DO 202 I=1,MX
      NJK(I)=0
202   XJK(I)=0.0
      MX=IMAX*KMAX
      DO 203 I=1,MX
      NIK(I)=0
203   XIK(I)=0.0
      DO 204 I=1,IMAX
      NI(I)=0
204   XI(I)=0.0
      DO 205 I=1,JMAX
      NJ(I)=0
205   XJ(I)=0.0
      DO 206 I=1,KMAX
      NK(I)=0
206   XK(I)=0.0
      SSTOT=0.0
      SUM=0.0
C
C     READ IN DATA
C
      IF(NWHICH.EQ.2)GO TO 5364
      IF(ISTD.EQ.0)GO TO 1
      FMT(1)='(3I,F'
      FMT(2)=')    '
1     READ(4,FMT,END=300)I,J,K,(X(L),L=1,NTIMES)
      IF(I.EQ.0)GO TO 300
      IF(I.GT.IMAX.OR.J.GT.JMAX.OR.K.GT.KMAX)GO TO 40
      IF(I.LT.1.OR.J.LT.1.OR.K.LT.1)GO TO 40
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
84    FORMAT(I)
      N=N+1
      NIJK(IJK)=NIJK(IJK)+1
      XIJK(IJK)=XIJK(IJK)+X(NTI)
	XXIJK(IJK)=XXIJK(IJK)+X(NTI)*X(NTI)
      SUM=SUM+X(NTI)
10    SSTOT=SSTOT+X(NTI)**2
      GO TO 1
5364  IF(ISTD.EQ.0)GO TO 3197
      FMT(1)='(10F)'
3197  READ(4,3198,END=300)(IENT(I),I=1,17)
3198  FORMAT(A1,15A5,A4)
      ENCODE(80,3198,IAST)(IENT(I),I=1,17)
      IF(IENT(1).EQ.'*')GO TO 3199
      IF(I.EQ.0)GO TO 3191
      DECODE(95,FMT,IAST)(X(L),L=1,NTIMES)
      N=N+1
      NIJK(IJK)=NIJK(IJK)+1
      XIJK(IJK)=XIJK(IJK)+X(NTI)
	XXIJK(IJK)=XXIJK(IJK)+X(NTI)*X(NTI)
      SUM=SUM+X(NTI)
      SSTOT=SSTOT+X(NTI)**2
      GO TO 3197
3199  DECODE(95,3140,IAST)I,J,K
3140   FORMAT(1X,3I)
      IF(I.EQ.0)GO TO 300
      IF(I.GT.IMAX.OR.J.GT.JMAX.OR.K.GT.KMAX)GO TO 3191
      IF(I.LT.1.OR.J.LT.1.OR.K.LT.1)GO TO 3191
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
      GO TO 3197
3191  TYPE 3192
3192  FORMAT( ' INDEX OUT OF RANGE USING METHOD 2. CANNOT CONTINUE!'/)
      CALL EXIT
40    TYPE 41
41    FORMAT(' INDEX OUT OF RANGE. OBSERVATION DELETED.'/)
      GO TO 1
300   CONTINUE
70    FORMAT('1',16A5)
      WRITE(6,70)(ID(I,NTI),I=1,5)
C
C      CALCULATE TOTAL AND WITHIN SS
C
      KIND=0
      MX=IMAX*JMAX*KMAX
      CT=SUM**2/N
C
      SSTOT=SSTOT-CT
C
      SUMM=0.0
      SSC=0.0
      NEM=0
      DO 91 I=1,MX
      IF(NIJK(I).EQ.0)GO TO 42
      SSC=SSC+XIJK(I)**2/NIJK(I)
      GO TO 91
42    NEM=NEM+1
91    CONTINUE
      SSC=SSC-CT
      SSABC=SSC
      SSW=SSTOT-SSC
      IF((N-MX).EQ.0)KIND=1
C
C          CALCULATE ACTUAL MARGINAL SUMS
C
      DO 900 I=1,IMAX
      DO 900 J=1,JMAX
      IJ=I+(J-1)*IMAX
      DO 900 K=1,KMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
C**** WMU-AM: #1.9.3, MOD=1, MTO, 19-OCT-77 ****
	IF(NIJK(IJK).GT.0)XXIJK(IJK)=XXIJK(IJK)-
	1				XIJK(IJK)*XIJK(IJK)/NIJK(IJK)
C**** END = MNL, #900-4 ****
      XIJ(IJ)=XIJ(IJ)+XIJK(IJK)
      NIJ(IJ)=NIJ(IJ)+NIJK(IJK)
      XI(I)=XI(I)+XIJK(IJK)
900   NI(I)=NI(I)+NIJK(IJK)
      DO 903 J=1,JMAX
      DO 903 K=1,KMAX
      JK=J+(K-1)*JMAX
      DO 903 I=1,IMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
	IF(NIJK(IJK).GT.1)XXIJK(IJK)=XXIJK(IJK)/(NIJK(IJK)-1)
      NJK(JK)=NJK(JK)+NIJK(IJK)
      XJK(JK)=XJK(JK)+XIJK(IJK)
      XJ(J)=XJ(J)+XIJK(IJK)
903   NJ(J)=NJ(J)+NIJK(IJK)
      DO 906 K=1,KMAX
      DO 906 I=1,IMAX
      IK=I+(K-1)*IMAX
      DO 906 J=1,JMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
      XIK(IK)=XIK(IK)+XIJK(IJK)
      NIK(IK)=NIK(IK)+NIJK(IJK)
      XK(K)=XK(K)+XIJK(IJK)
 906  NK(K)=NK(K)+NIJK(IJK)
      IF(IFM.EQ.'YES')CALL PRINT(IMAX,JMAX,KMAX,XI,XJ,XK,XIJ,XJK,XIK,
     1XIJK,NI,NJ,NK,NIJ,NJK,NIK,NIJK,SUM,N,INAME,JNAME,KNAME,XXIJK)
	CALL BARTLT(IMAX,JMAX,KMAX,NIJK,XXIJK,BART)
	IDEG=IMAX*JMAX*KMAX-1
	IF(BART.LT.0)GO TO 9127
C---------------CBART, IDEG INPUT AND YPROB RETURNED.
	CALL CHIPRB(BART,IDEG,YPROB)
	IF(IERR.NE.99)WRITE(6,9128)YPROB
9128	FORMAT(5X,'WITH CHI-SQUARE PROBABILITY=',T39,F5.3)
	IF(IERR.EQ.99)WRITE(6,9129)
9129	FORMAT(5X,'NO PROBABILITY CALC. DF OR BARTLETT OUTSIDE LIMITS')
9127      KON=NIJK(1)
      DO 1012 I=2,IMAX*JMAX*KMAX
      IF(NIJK(I).NE.KON)GO TO 1014
1012  CONTINUE
      GO TO 1020
C**** WMU-AM: #1.9.3, MOD=1, MTO, 19-SEP-77 ****
1014	DO 2000 I=1,IMAX
	DO 2000 J=1,JMAX
	DO 2000 K=1,KMAX
	IJK=I+(J-1+(K-1)*JMAX)*IMAX
	IF(NIJK(IJK).NE.NI(I)*NJ(J)*NK(K)/(N*N))GOTO 1018
2000	CONTINUE
	CALL PROP(IMAX,JMAX,KMAX,XI,XJ,XK,XIJ,XIK,XJK,XIJK,
	1	NI,NJ,NK,NIJ,NJK,NIK,NIJK,
	2	SSTOT,SSW,N,CT,SSC)
1018  CALL FIT(IMAX,JMAX,KMAX,XI,XJ,XK,NI,NJ,NK,XIJ,XJK,XIK,
     1NIJ,NJK,NIK,NIJK,D,DF,XD,FF,YB,A,B,C,MC,MR,SUM,SSW,N,SSTOT,NEM,
     2SSABC,INAME,JNAME,KNAME,SSINT,CT)
C**** END = MNL, #1016-2 ****
      DO 1016 I=1,IMAX*JMAX*KMAX
      IF(NIJK(I).EQ.0)GO TO 1017
1016  CONTINUE
      CALL WEIGH(IMAX,JMAX,KMAX,NIJK,XIJK,INAME,JNAME,KNAME,A,B,C,YB,XD,
     1DF,SSINT,SSTOT,SSW,N)
      WRITE(6,1021)
1021  FORMAT(///,20X,'UNWEIGHTED MEANS ANOVA')
      GO TO 1023
1020  WRITE(6,1022)
1022  FORMAT(///,20X,'ANALYSIS OF VARIANCE')
1023  CALL UNWEI(IMAX,JMAX,KMAX,XI,XJ,XK,XIJ,XJK,XIK,XIJK,NI,NJ,NK,NIJ,
     1NJK,NIK,NIJK,SSTOT,SSW,N,INAME,JNAME,KNAME)
1017  IF(NTI.NE.NTIMES)GO TO 81
      RETURN
81    REWIND 5
      GO TO 80
      END
C
C		BARTLETTS TEST FOR HOMOGENEITY
C
C---------------IMAX, JMAX, KMAX, NIJK, XXIJK ARE INPUT;  CS RETURNED
	SUBROUTINE BARTLT(IMAX,JMAX,KMAX,NIJK,XXIJK,CS)
	DIMENSION NIJK(1),XXIJK(1)
	V=0
	VI=0
	SPSQ=0
	CS=0
	WRITE(6,3)
3	FORMAT(///,2X,'BARTLETT''S TEST STATISTIC FOR TESTING ',
     1	'HOMOGENEITY OF CELL VARIANCES')
	DO 1 I=1,IMAX
	DO 1 J=1,JMAX
	DO 1 K=1,KMAX
	IJK=I+(J-1+(K-1)*JMAX)*IMAX
	IF(NIJK(IJK).LT.2)GO TO 99
	V=V+(NIJK(IJK)-1)
	VI=VI+1./(NIJK(IJK)-1)
	CS=CS+(NIJK(IJK)-1)*ALOG(XXIJK(IJK))
1	SPSQ=SPSQ+XXIJK(IJK)*(NIJK(IJK)-1)
	SPSQ=SPSQ/V
	IJKMAX=IMAX*JMAX*KMAX
	ID=IJKMAX-1
	CS=(V*ALOG(SPSQ)-CS)/(1.+(VI-1./V)/(3.*ID))
	WRITE(6,4)IJKMAX,CS,ID
4	FORMAT(/5X,'NUMBER OF CELL VARIANCES='I4,/,
     1	5X,'BARTLETT''S STATISTIC=',F12.3,5X,'DF=',I4)
	RETURN
99	WRITE(6,5)
5	FORMAT('  IS NOT POSSIBLE SINCE AT LEAST ONE SAMPLE SIZE',
     1	' IS EITHER 0 OR 1')
	CS=-1
	RETURN
	END
C
C		CUMULATIVE CHI SQUARE SUBROUTINE
C
C		SOURCE-CHARLES A. NAGY
C
C		CUMULATIVE CHI SQUARE SUBROUTINE
C
C		SOURCE-CHARLES A. NAGY
C
C          FITTING CONSTANTS
C
C**** WMU-AM: #1.9.3, MOD=1, MTO, 23-SEP-77 ****
C---------------IMAX, JMAX, KMAX, XI, XJ, XK, NI, NJ, NK, XIJ, XJK,
C---------------XIK, NIJ, NJK, NIK, NIJK, INAME, JNAME, KNAME, SSTOT, SSW,
C--------------- N, CT, NEM, SSABC ARE INPUT.  D, DF, XD, FF, YB,
C--------------- A, B, C, SSINT ARE RETURNED.  SPACE FOR MC, MR IS RESERVED
C--------------- BY ALLCOR AND USED BY MINVSQ BUT NOT NEEDED ELSEWHERE
      SUBROUTINE FIT(IMAX,JMAX,KMAX,XI,XJ,XK,NI,NJ,NK,XIJ,XJK,XIK,
     1NIJ,NJK,NIK,NIJK,D,DF,XD,FF,YB,A,B,C,MC,MR,SUM,SSW,N,SSTOT,NEM,
     2SSABC,INAME, JNAME,KNAME,SSINT,CT)
      DIMENSION XI(1),XJ(1),XK(1),NI(1),NJ(1),NK(1),XIJ(1),XJK(1),XIK(1)
      DIMENSION NIJ(1),NJK(1),NIK(1),NIJK(1),D(1),DF(1),FF(1),YB(1)
      DIMENSION XD(1),A(1),B(1),C(1),MC(1),MR(1)
C**** END = FIT, #1104-10 ****
      JM1=JMAX+1
      DO 1100 I=1,JM1
      DO 1100 J=1,JM1
      IJ=I+(J-1)*JM1
      IF(I.EQ.JM1.OR.J.EQ.JM1)GO TO 1103
      XNTL=0.0
      DO 1104 K=1,IMAX
      LM1=K+(I-1)*IMAX
      LM2=K+(J-1)*IMAX
      XNTL=XNTL+FLOAT(NIJ(LM1)*NIJ(LM2))/FLOAT(NI(K))
1104   CONTINUE
      IF(I.EQ.J)GO TO 1102
      D(IJ)=-XNTL
      GO TO 1100
1102   D(IJ)=NJ(J)-XNTL
      GO TO 1100
1103  D(IJ)=1.0
      IF(I.EQ.J)D(IJ)=0.0
1100  CONTINUE
      CALL MINVSQ(D,JM1,1.0,MC,MR,JM1,-1,2,DET,IEXP)
      DO 1105 J=1,JMAX
      XNTL=0.0
      DO 1106 I=1,IMAX
      IJ=I+(J-1)*IMAX
1106  XNTL=XNTL+FLOAT(NIJ(IJ))*XI(I)/FLOAT(NI(I))
1105  YB(J)=XNTL
      DO 1107 J=1,JMAX
      XNTL=0.0
      DO 1108 K=1,JMAX
      JK=J+(K-1)*JM1
1108  XNTL=XNTL+D(JK)*(XJ(K)-YB(K))
1107   XD(J)=XNTL
      DO 1109 J=1,JMAX
      DO 1110 K=1,KMAX
      XNTL=0.0
      JK=J+(K-1)*JMAX
      DO 1111I=1,IMAX
      IJ=I+(J-1)*IMAX
      IK=I+(K-1)*IMAX
1111  XNTL=XNTL+FLOAT(NIJ(IJ)*NIK(IK))/FLOAT(NI(I))
1110  DF(JK)=FLOAT(NJK(JK))-XNTL
1109  CONTINUE
      DO 1112 J=1,JMAX
      DO 1113 K=1,KMAX
      XNTL=0.0
      JK=J+(K-1)*JMAX
      DO 1114 L=1,JMAX
      JL=J+(L-1)*JM1
      LK=L+(K-1)*JMAX
1114  XNTL=XNTL+D(JL)*DF(LK)
1113  FF(JK)=XNTL
1112  CONTINUE
      DO 1115 K=1,KMAX
      XNTL=0.0
      DO 1116 I=1,IMAX
      IK=I+(K-1)*IMAX
1116  XNTL=XNTL+FLOAT(NIK(IK))*XI(I)/FLOAT(NI(I))
      YNTL=0.0
      DO 1117 J=1,JMAX
      JK=J+(K-1)*JMAX
1117  YNTL=YNTL+DF(JK)*XD(J)
1115  YB(K)=XK(K)-XNTL-YNTL
      KM1=KMAX+1
      DO 1118 K=1,KM1
      DO 1118 KK=1,KM1
      KKK=K+(KK-1)*KM1
      IF(K.EQ.KM1)GO TO 1119
      IF(KK.EQ.KM1)GO TO 1119
      XNTL=0.0
      DO 1120 I=1,IMAX
      IK=I+(K-1)*IMAX
      IKK=I+(KK-1)*IMAX
1120  XNTL=XNTL+FLOAT(NIK(IK)*NIK(IKK))/FLOAT(NI(I))
      YNTL=0.0
      DO 1122  J=1,JMAX
      JK=J+(K-1)*JMAX
      JKK=J+(KK-1)*JMAX
1122  YNTL=YNTL+DF(JK)*FF(JKK)
      IF(K.EQ.KK)GO TO 1121
      D(KKK)=-XNTL-YNTL
      GO TO 1118
1121  D(KKK)=NK(K)-XNTL-YNTL
      GO TO 1118
1119  D(KKK)=1.0
      IF(K.EQ.KK)D(KKK)=0.0
1118  CONTINUE
      CALL MINVSQ(D,KM1,1.0,MC,MR,KM1,-1,2,DET,IEXP)
      CTERM=0.0
      DO 1123 K=1,KMAX
      XNTL=0.0
      DO 1124 KK=1,KMAX
      KKK=K+(KK-1)*KM1
1124  XNTL=XNTL+D(KKK)*YB(KK)
      C(K)=XNTL
1123  CTERM=CTERM+C(K)*XK(K)
      BTERM=0.0
      DO 1125 J=1,JMAX
      XNTL=0.0
      DO 1126 K=1,KMAX
      JK=J+(K-1)*JMAX
1126  XNTL=XNTL+FF(JK)*C(K)
      B(J)=XD(J)-XNTL
1125  BTERM=BTERM+B(J)*XJ(J)
      ZNTL=0.0
      DO 1127 I=1,IMAX
      XNTL=0.0
      DO 1128 J=1,JMAX
      IJ=I+(J-1)*IMAX
1128  XNTL=XNTL+FLOAT(NIJ(IJ))*B(J)
      YNTL=0.0
      DO 1129 K=1,KMAX
      IK=I+(K-1)*IMAX
1129  YNTL=YNTL+FLOAT(NIK(IK))*C(K)
      A(I)=(XI(I)-XNTL-YNTL)/FLOAT(NI(I))
1127  ZNTL=ZNTL+A(I)
      XMU=ZNTL/FLOAT(IMAX)
      ATERM=0.0
      DO 1130 I=1,IMAX
1130   ATERM=ATERM+(A(I)-XMU)*XI(I)
      RABC=XMU*SUM+ATERM+BTERM+CTERM-CT
      SSINT=SSABC-RABC
      ZNTL=0.0
      DO 1140 I=1,IMAX
      XNTL=0.0
      DO 1141 J=1,JMAX
      IJ=I+(J-1)*IMAX
      XNTL=XNTL+FLOAT(NIJ(IJ))*XD(J)/FLOAT(NI(I))
1141  CONTINUE
      YB(I)=XI(I)/FLOAT(NI(I))-XNTL
      ZNTL=ZNTL+YB(I)
1140  CONTINUE
      ZNTL=ZNTL/FLOAT(IMAX)
      RAB=ZNTL*SUM-CT
      DO 1142 I=1,IMAX
      RAB=RAB+(YB(I)-ZNTL)*XI(I)
1142  CONTINUE
      DO 1143 J=1,JMAX
      RAB=RAB+XD(J)*XJ(J)
1143  CONTINUE
      RC=RABC-RAB
      KM1=KMAX+1
      DO 1144 K=1,KM1
      DO 1144 L=1,KM1
      KL=K+(L-1)*KM1
      IF(K.EQ.L)GO TO 1145
      IF(K.EQ.KM1)GO TO 1146
      IF(L.EQ.KM1)GO TO 1146
      XNTL=0.0
      DO 1147 I=1,IMAX
      IK=I+(K-1)*IMAX
      IL=I+(L-1)*IMAX
      XNTL=XNTL+FLOAT(NIK(IK)*NIK(IL))/FLOAT(NI(I))
1147  CONTINUE
      D(KL)=-XNTL
      GO TO 1144
1145  IF(K.EQ.KM1)GO TO 1245
      XNTL=0.0
      DO 1244 I=1,IMAX
      IK=I+(K-1)*IMAX
      XNTL=XNTL+FLOAT(NIK(IK))**2/FLOAT(NI(I))
1244  CONTINUE
      D(KL)=FLOAT(NK(K))-XNTL
      GO TO 1144
1245  D(KL)=0.0
      GO TO 1144
1146  D(KL)=1.0
1144  CONTINUE
      CALL MINVSQ(D,KM1,1.0,MC,MR,KM1,-1,2,DET,IEXP)
      DO 1148 K=1,KMAX
      XNTL=0.0
      DO 1149 I=1,IMAX
      IK=I+(K-1)*IMAX
      XNTL=XNTL+FLOAT(NIK(IK))*XI(I)/FLOAT(NI(I))
1149  CONTINUE
      YB(K)=XNTL
1148  CONTINUE
      DO 1150 K=1,KMAX
      XNTL=0.0
      DO 1151 KK=1,KMAX
      KKK=K+(KK-1)*KM1
      XNTL=XNTL+D(KKK)*(XK(KK)-YB(KK))
1151  CONTINUE
      XD(K)=XNTL
1150  CONTINUE
      ZNTL=0.0
      DO 1152 I=1,IMAX
      XNTL=0.0
      DO 1153 K=1,KMAX
      IK=I+(K-1)*IMAX
      XNTL=XNTL+FLOAT(NIK(IK))*XD(K)/FLOAT(NI(I))
1153  CONTINUE
      YB(I)=XI(I)/FLOAT(NI(I))-XNTL
      ZNTL=ZNTL+YB(I)
1152  CONTINUE
      ZNTL=ZNTL/FLOAT(IMAX)
      RAC=ZNTL*SUM-CT
      DO 1154 I=1,IMAX
      RAC=RAC+(YB(I)-ZNTL)*XI(I)
1154  CONTINUE
      DO 1155 K=1,KMAX
      RAC=RAC+XD(K)*XK(K)
1155  CONTINUE
      RB=RABC-RAC
      DO 1156 K=1,KM1
      DO 1157 L=1,KM1
      KL=K+(L-1)*KM1
      IF(K.EQ.L)GO TO 1158
      IF(K.EQ.KM1)GO TO 1159
      IF(L.EQ.KM1)GO TO 1159
      XNTL=0.0
      DO 1160 J=1,JMAX
      JK=J+(K-1)*JMAX
      JL=J+(L-1)*JMAX
      XNTL=XNTL+FLOAT(NJK(JK)*NJK(JL))/FLOAT(NJ(J))
1160  CONTINUE
      D(KL)=-XNTL
      GO TO 1157
1158   IF(K.EQ.KM1)GO TO 1161
      XNTL=0.0
      DO 1162 J=1,JMAX
      JK=J+(K-1)*JMAX
      XNTL=XNTL+FLOAT(NJK(JK))**2/FLOAT(NJ(J))
1162  CONTINUE
      D(KL)=FLOAT(NK(K))-XNTL
      GO TO 1157
1161  D(KL)=0.0
      GO TO 1157
1159  D(KL)=1.0
1157  CONTINUE
1156  CONTINUE
      CALL MINVSQ(D,KM1,1.0,MC,MR,KM1,-1,2,DET,IEXP)
      DO 1163 K=1,KMAX
      XNTL=0.0
      DO 1164 J=1,JMAX
      JK=J+(K-1)*JMAX
      XNTL=XNTL+FLOAT(NJK(JK))*XJ(J)/FLOAT(NJ(J))
1164  CONTINUE
      YB(K)=XNTL
1163  CONTINUE
      DO 1165 K=1,KMAX
      XNTL=0.0
      DO 1166 L=1,KMAX
      KL=K+(L-1)*KM1
      XNTL=XNTL+D(KL)*(XK(L)-YB(L))
1166  CONTINUE
      XD(K)=XNTL
1165  CONTINUE
      ZNTL=0.0
      DO 1167 J=1,JMAX
      XNTL=0.0
      DO 1168 K=1,KMAX
      JK=J+(K-1)*JMAX
      XNTL=XNTL+FLOAT(NJK(JK))*XD(K)/FLOAT(NJ(J))
 1168 CONTINUE
      YB(J)=XJ(J)/FLOAT(NJ(J))-XNTL
      ZNTL=ZNTL+YB(J)
1167  CONTINUE
      ZNTL=ZNTL/FLOAT(JMAX)
      RBC=ZNTL*SUM-CT
      DO 1171 J=1,JMAX
      RBC=RBC+(YB(J)-ZNTL)*XJ(J)
1171  CONTINUE
      DO 1169 K=1,KMAX
      RBC=RBC+XD(K)*XK(K)
1169  CONTINUE
      RA=RABC-RBC
      WRITE(6,100)
100   FORMAT(///,17X,'PRELIMINARY ANOVA - FITTING CONSTANTS')
      WRITE(6,101)
101   FORMAT(/5X,'SOURCE',5X,'SUM OF SQUARES',3X,'DF',4X,'MEAN SQUARE',
     14X,'F',4X,'PROB')
      IDF=IMAX*JMAX*KMAX-1
C**** WMU-AM: #1.9.3, MOD=1, MTO, 19-SEP-77 ****
	DO 201 I=1,IMAX
	DO 201 J=1,JMAX
	DO 201 K=1,KMAX
	IJK=I+(J-1+(K-1)*JMAX)*IMAX
201	IF(NIJK(IJK).EQ.0) IDF=IDF-1
C**** WMU-AM: #1.9.3, MOD=1, MTO, 7-OCT-77 ****
	XMS=SSABC/FLOAT(IDF)
	KDF=(IMAX*JMAX*KMAX-1)-IDF
C**** END = FIT, #202-6 ****
	IDW=N-IDF-1
	IF(IDF.GT.0.AND.IDW.GT.0) GOTO 202
	TYPE 2001
2001	FORMAT (' NOT ENOUGH DATA TO CARRY OUT AN ANOVA.')
	WRITE (6,2001)
	RETURN
202   SMW=SSW/FLOAT(IDW)
C**** END = FIT, #102-3 ****
      F=XMS/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,102)SSABC,IDF,XMS,F,PROB
102   FORMAT(/5X,'SUBCLASS',2X,F15.2,I5,F15.2,F7.2,F6.3)
      IDF=IMAX+JMAX+KMAX-3
      WRITE(6,103)RABC,IDF
103   FORMAT(/3X,'MAIN EFFECTS',F15.2,I5)
      IDF=IMAX*JMAX*KMAX-IMAX-JMAX-KMAX+2-NEM
C**** WMU-AM: #1.9.3, MOD=1, MTO, 10-OCT-77 ****
	IF(IDF.GT.0) GOTO 204
	WRITE(6,2002)
2002	FORMAT (' TOO MANY MISSING CELLS TO OBTAIN AN INTER-'/
	1	' ACTION TEST.'//)
	GOTO 205
204   XMS=SSINT/FLOAT(IDF)
      F=XMS/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,104)SSINT,IDF,XMS,F,PROB
104   FORMAT(/4X'INTERACTION',F15.2,I5,F15.2,F7.2,F6.3)
205   WRITE(6,105)SSW,IDW,SMW
C**** END = FIT, #105-1 ****
105   FORMAT(/6X,'WITHIN',3X,F15.2,I5,F15.2)
      IDF=N-1
      WRITE(6,106)SSTOT,IDF
106   FORMAT(/6X,'TOTAL',4X,F15.2,I5)
      WRITE(6,107)
107   FORMAT(///,20X,'FINAL ANOVA - FITTING CONSTANTS')
      WRITE(6,108)
108   FORMAT(/7X,'SOURCE',8X,'SUM OF SQUARES    DF   MEAN SQUARE     F'
     13X,'PROB')
      IDT=IMAX+JMAX+KMAX-3
      WRITE(6,109)RABC,IDT
109   FORMAT(//4X,'MAIN EFFECTS',4X,F15.2,I5)
      IDF=IDT-(IMAX-1)
      WRITE(6,110)JNAME,KNAME,INAME,RBC,IDF
110   FORMAT(//1X,A1,' & ',A1,' IGNORING ',A1,3X,F15.2,I5)
      IDF=IMAX-1
      XMS=RA/FLOAT(IDF)
      F=XMS/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,111)INAME,JNAME,KNAME,RA,IDF,XMS,F,PROB
111   FORMAT(/1X,A1,' ELIMINATING ',A1,' & ',A1,F15.2,
	1	I5,F15.2,F7.2,F6.3)
      IDF=IDT-(JMAX-1)
      WRITE(6,110)INAME,KNAME,JNAME,RAC,IDF
      IDF=JMAX-1
      XMS=RB/FLOAT(IDF)
      F=XMS/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,111)JNAME,INAME,KNAME,RB,IDF,XMS,F,PROB
      IDF=IDT-(KMAX-1)
      WRITE(6,110)INAME,JNAME,KNAME,RAB,IDF
      IDF=KMAX-1
      XMS=RC/FLOAT(IDF)
      F=XMS/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,111)KNAME,INAME,JNAME,RC,IDF,XMS,F,PROB
      WRITE(6,112)SSW,IDW,SMW
112   FORMAT(//7X,'WITHIN',7X,F15.2,I5,F15.2)
      RETURN
      END
C
C            PRINT OUT MEANS
C
C---------------ALL ARGUMENTS ARE INPUT.
      SUBROUTINE PRINT(IMAX,JMAX,KMAX,XI,XJ,XK,XIJ,XJK,XIK,XIJK,NI,
     1NJ,NK,NIJ,NJK,NIK,NIJK,SUM,N,INAME,JNAME,KNAME,XXIJK)
      DIMENSION XI(1),XJ(1),XK(1),XIJ(1),XJK(1),XIK(1),XIJK(1),NI(1)
     1,NJ(1),NK(1),NIJ(1),NJK(1),NIK(1),NIJK(1),XXIJK(1)
311   GMEAN=SUM/FLOAT(N)
      WRITE(6,601)INAME,JNAME ,KNAME
601   FORMAT(///,10X,A1,3X,A1,3X,A1,6X,'MEAN',13X,'STD DEV',5X,'N'/)
      DO 602 I=1,IMAX
      DO 602 J=1,JMAX
      DO 602 K=1,KMAX
      IND=I+(J-1+(K-1)*JMAX)*IMAX
      IF(NIJK(IND).EQ.0)GO TO 667
	STDV=SQRT(XXIJK(IND))
      XX=XIJK(IND)/FLOAT(NIJK(IND))
      WRITE(6,603)I,J,K,XX,STDV,NIJK(IND)
603   FORMAT(9X,I3,1X,I3,1X,I3,G,F15.7,I5)
      GO TO 602
667   WRITE(6,668)I,J,K,NIJK(IND)
668   FORMAT(9X,I3,1X,I3,1X,I3,30X,I5)
602   CONTINUE
      DO 604 I=1,IMAX
      DO 604 J=1,JMAX
      IND=I+(J-1)*IMAX
      XX=XIJ(IND)/FLOAT(NIJ(IND))
604   WRITE(6,605)I,J,XX,NIJ(IND)
605   FORMAT(9X,I3,1X,I3,1X,'  .',G,15X,I5)
      DO 606 I=1,IMAX
      DO 606 K=1,KMAX
      IND=I+(K-1)*IMAX
      XX=XIK(IND)/FLOAT(NIK(IND))
606   WRITE(6,607)I,K,XX,NIK(IND)
607   FORMAT(9X,I3,1X,'  .',1X,I3,G,15X,I5)
      DO 608 J=1,JMAX
      DO 608 K=1,KMAX
      IND=J+(K-1)*JMAX
      XX=XJK(IND)/FLOAT(NJK(IND))
608   WRITE(6,609)J,K,XX,NJK(IND)
609   FORMAT(9X,'  .',1X,I3,1X,I3,G,15X,I5)
      DO 610 I=1,IMAX
      XX=XI(I)/FLOAT(NI(I))
610   WRITE(6,611)I,XX,NI(I)
611   FORMAT(9X,I3,1X,'  .',1X,'  .',G,15X,I5)
      DO 612 J=1,JMAX
      XX=XJ(J)/FLOAT(NJ(J))
612   WRITE(6,613)J,XX,NJ(J)
613   FORMAT(9X,'  .',1X,I3,1X,'  .',G,15X,I5)
      DO 614 K=1,KMAX
      XX=XK(K)/FLOAT(NK(K))
614   WRITE(6,615)K,XX,NK(K)
615   FORMAT(9X,'  .',1X,'  .',1X,I3,G,15X,I5)
      WRITE(6,616)GMEAN,N
616   FORMAT(9X,'  .',1X,'  .',1X,'  .',G,15X,I5)
      RETURN
      END
C
C           WEIGHTED MEANS ANALYSIS
C
C---------------YS, XN1, XN2, XN3, YT, YR ARE RETURNED OTHER ARGS.
C--------------- ARE INPUT.
      SUBROUTINE WEIGH(IMAX,JMAX,KMAX,NIJK,XIJK,INAME,JNAME,KNAME,
     1XN1,XN2,XN3,YR,YS,YT,SSINT,SSTOT,SE,N)
      DIMENSION NIJK(1),XIJK(1),XN1(1),XN2(1),XN3(1),YR(1),YS(1),YT(1)
      XN=0.0
      XD=0.0
      DO 1 I=1,IMAX
      XNTL=0.0
      ZNTL=0.0
      DO 2 J=1,JMAX
      DO 3 K=1,KMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
      XNTL=XNTL+1.0/FLOAT(NIJK(IJK))
      ZNTL=ZNTL+XIJK(IJK)/NIJK(IJK)
3     CONTINUE
2     CONTINUE
      XN1(I)=FLOAT(JMAX*JMAX*KMAX*KMAX)/XNTL
      YR(I)=ZNTL/FLOAT(JMAX*KMAX)
      XN=XN+XN1(I)*YR(I)
      XD=XD+XN1(I)
1     CONTINUE
      CCR=XN/XD
      XN=0.0
      XD=0.0
      DO 4 J=1,JMAX
      ZNTL=0.0
      XNTL=0.0
      DO 5 I=1,IMAX
      DO 6 K=1,KMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
      XNTL=XNTL+1.0/FLOAT(NIJK(IJK))
      ZNTL=ZNTL+XIJK(IJK)/FLOAT(NIJK(IJK))
 6    CONTINUE
5     CONTINUE
      XN2(J)=FLOAT(IMAX*IMAX*KMAX*KMAX)/XNTL
      YS(J)=ZNTL/FLOAT(IMAX*KMAX)
      XN=XN+XN2(J)*YS(J)
      XD=XD+XN2(J)
4     CONTINUE
      CCS=XN/XD
      XN=0.0
      XD=0.0
      DO 7 K=1,KMAX
      ZNTL=0.0
      XNTL=0.0
      DO 8 I=1,IMAX
      DO 9 J=1,JMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
      XNTL=XNTL+1.0/FLOAT(NIJK(IJK))
      ZNTL=ZNTL+XIJK(IJK)/FLOAT(NIJK(IJK))
9     CONTINUE
8      CONTINUE
      XN3(K)=FLOAT(IMAX*IMAX*JMAX*JMAX)/XNTL
      YT(K)=ZNTL/FLOAT(IMAX*JMAX)
      XN=XN+XN3(K)*YT(K)
      XD=XD+XN3(K)
7     CONTINUE
      CCT=XN/XD
      SSR=0.0
      DO 10 I=1,IMAX
10    SSR=SSR+XN1(I)*((YR(I)-CCR)**2)
      SSS=0.0
      DO 11 J=1,JMAX
11    SSS=SSS+XN2(J)*((YS(J)-CCS)**2)
      SST=0.0
      DO 12 K=1,KMAX
12    SST=SST+XN3(K)*((YT(K)-CCT)**2)
      WRITE(6,100)
100   FORMAT(///,22X,'THREE WAY WEIGHTED MEANS ANOVA')
      WRITE(6,101)
101   FORMAT(/4X,'SOURCE',9X,'SUM OF SQUARES  DF',4X,'MEAN SQUARE',5X,
     1'F   PROB')
      IDF=IMAX-1
      XMS=SSR/FLOAT(IDF)
      IDW=N-IMAX*JMAX*KMAX
      SMW=SE/FLOAT(IDW)
      F=XMS/SMW
      PROB=FISHER (IDF,IDW,F)
      WRITE(6,102)INAME,SSR,IDF,XMS,F,PROB
102   FORMAT(/6X,A1,10X,F15.2,I5,F15.2,F7.2,F6.3)
      IDF=JMAX-1
      XMS=SSS/FLOAT(IDF)
      F=XMS/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,102)JNAME,SSS,IDF,XMS,F,PROB
      IDF=KMAX-1
      XMS=SST/FLOAT(IDF)
      F=XMS/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,102)KNAME,SST,IDF,XMS,F,PROB
      IDF=IMAX*JMAX*KMAX-IMAX-JMAX-KMAX+2
      XMS=SSINT/FLOAT(IDF)
      F=XMS/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,103)SSINT,IDF,XMS,F,PROB
103   FORMAT(/3X,'INTERACTION',3X,F15.2,I5,F15.2,F7.2,F6.3)
      WRITE(6,104)SE,IDW,SMW
104   FORMAT(/5X,'WITHIN',6X,F15.2,I5,F15.2)
      IDF=N-1
      WRITE(6,105)SSTOT,IDF
105   FORMAT(/5X,'TOTAL',7X,F15.2,I5)
      RETURN
      END
C
C     UNWEIGHTED MEANS ANALYSIS
C
C---------------XI, XJ, XK, XIJ, XJK, XIK ARE MODIFIED.  OTHER
C--------------- ARGS. ARE INPUT WITHOUT MODIFICATION.
      SUBROUTINE UNWEI(IMAX,JMAX,KMAX,XI,XJ,XK,XIJ,XJK,XIK,XIJK,
     1NI,NJ,NK,NIJ,NJK,NIK,NIJK,SSTOT,SSW,N,INAME,JNAME,KNAME)
      DIMENSION XI(1),XJ(1),XK(1),XIJ(1),XJK(1),XIK(1),XIJK(1),NI(1),
     1NJ(1),NK(1),NIJ(1),NJK(1),NIK(1),NIJK(1)
      MX=IMAX*JMAX*KMAX
      XN=0.0
      DO 100 I=1,MX
      XN=XN+1.0/FLOAT(NIJK(I))
100   XIJK(I)=XIJK(I)/FLOAT(NIJK(I))
      XN=XN/FLOAT(MX)
      KIND=0
      IF(N.EQ.MX)KIND=1
      SUM=0.0
      DO 309 I=1,MX
309   SUM=SUM+XIJK(I)
      CTT=SUM**2/MX
C
C      REINITIALIZE CERTAIN VALUES
C
      DO 909 I=1,IMAX
909   XI(I)=0.0
      DO 910 J=1,JMAX
910   XJ(J)=0.0
      DO 911 K=1,KMAX
911   XK(K)=0.0
      DO 912 I=1,IMAX*JMAX
912   XIJ(I)=0.0
      DO 913 I=1,JMAX*KMAX
913   XJK(I)=0.0
      DO 914 I=1,IMAX*KMAX
914   XIK(I)=0.0
      SSI=0.0
      SSJ=0.0
      SSK=0.0
      SSIJ=0.0
      SSIK=0.0
      SSJK=0.0
      SSIJK=0.0
C
C     CALCULATE MARGINAL SUMS
C
      DO 301 I=1,IMAX
      DO 301 J=1,JMAX
      IJ=I+(J-1)*IMAX
      DO 302 K=1,KMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
      XIJ(IJ)=XIJ(IJ)+XIJK(IJK)
302   XI(I)=XI(I)+XIJK(IJK)
301   CONTINUE
      DO 303 J=1,JMAX
      DO 304 K=1,KMAX
      JK=J+(K-1)*JMAX
      DO 305 I=1,IMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
      XJK(JK)=XJK(JK)+XIJK(IJK)
305   XJ(J)=XJ(J)+XIJK(IJK)
304   CONTINUE
303   CONTINUE
      DO 306 K=1,KMAX
      DO 307 I=1,IMAX
      IK=I+(K-1)*IMAX
      DO 308 J=1,JMAX
      IJK=I+(J-1+(K-1)*JMAX)*IMAX
      XIK(IK)=XIK(IK)+XIJK(IJK)
308   XK(K)=XK(K)+XIJK(IJK)
307   CONTINUE
306   CONTINUE
C
C     COMPUTE MAIN EFFECTS AND INTERACTIONS
C
      SSC=0.0
      DO 500 I=1,MX
500   SSC=SSC+XIJK(I)**2
      MX=IMAX*JMAX
      DO 501 I=1,MX
501   SSIJ=SSIJ+XIJ(I)**2
      MX=JMAX*KMAX
      DO 502 I=1,MX
502   SSJK=SSJK+XJK(I)**2
      MX=IMAX*KMAX
      DO 503 I=1,MX
503   SSIK=SSIK+XIK(I)**2
      DO 504 I=1,IMAX
504   SSI=SSI+XI(I)**2
      DO 505 I=1,JMAX
505   SSJ=SSJ+XJ(I)**2
      DO 506 I=1,KMAX
506   SSK=SSK+XK(I)**2
      XN=1.0/XN
      SSC=(SSC-CTT)*XN
      SSI=(SSI/FLOAT(JMAX*KMAX)-CTT)*XN
      SSJ=(SSJ/FLOAT(IMAX*KMAX)-CTT)*XN
      SSK=(SSK/FLOAT(IMAX*JMAX)-CTT)*XN
      SSIJ=(SSIJ/FLOAT(KMAX)-CTT)*XN-SSI-SSJ
      SSJK=(SSJK/FLOAT(IMAX)-CTT)*XN-SSJ-SSK
      SSIK=(SSIK/FLOAT(JMAX)-CTT)*XN-SSI-SSK
      SSIJK=SSC-SSI-SSJ-SSK-SSIJ-SSIK-SSJK
C
C      CALCULATE AND PRINT ANOVA
C
      WRITE(6,700)XN
700   FORMAT(/20X,'H.M. CELL SIZE = ',F10.5)
      IDW=N-IMAX*JMAX*KMAX
      IF(KIND.EQ.1)IDW=(IMAX-1)*(JMAX-1)*(KMAX-1)
      SS=SSW
      IF(KIND.EQ.1)SS=SSIJK
      SMW=SS/FLOAT(IDW)
      WRITE(6,702)
702   FORMAT(/7X,'SOURCE',4X,'SUM OF SQUARES',4X,'DF',3X,'MEAN SQUARE',7
     1X'F',5X,'PROB')
      IDF=IMAX*JMAX*KMAX-1
      SM=SSC/FLOAT(IDF)
      F=SM/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,703)SSC,IDF,SM,F,PROB
703   FORMAT(/4X,'SUBCLASSES',2X,F15.2,I6,F15.2,F9.3,F8.3)
      IDF=IMAX-1
      SM=SSI/FLOAT(IDF)
      F=SM/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,704)INAME,SSI,IDF,SM,F,PROB
 704  FORMAT(/8X,A1,7X,F15.2,I6,F15.2,F9.3,F8.3)
      IDF=JMAX-1
      SM=SSJ/FLOAT(IDF)
      F=SM/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,704)JNAME,SSJ,IDF,SM,F,PROB
      IDF=KMAX-1
      SM=SSK/FLOAT(IDF)
      F=SM/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,704)KNAME,SSK,IDF,SM,F,PROB
      IDF=(IMAX-1)*(JMAX-1)
      SM=SSIJ/FLOAT(IDF)
      F=SM/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,705)INAME,JNAME,SSIJ,IDF,SM,F,PROB
705   FORMAT(/8X,2A1,6X,F15.2,I6,F15.2,F9.3,F8.3)
      IDF=(IMAX-1)*(KMAX-1)
      SM=SSIK/FLOAT(IDF)
      F=SM/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,705)INAME,KNAME,SSIK,IDF,SM,F,PROB
      IDF=(JMAX-1)*(KMAX-1)
      SM=SSJK/FLOAT(IDF)
      F=SM/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,705)JNAME,KNAME,SSJK,IDF,SM,F,PROB
      IF(KIND.EQ.1)GO TO 710
      IDF=(IMAX-1)*(JMAX-1)*(KMAX-1)
      SM=SSIJK/FLOAT(IDF)
      F=SM/SMW
      PROB=FISHER(IDF,IDW,F)
      WRITE(6,706)INAME,JNAME,KNAME,SSIJK,IDF,SM,F,PROB
706    FORMAT(/8X,3A1,5X,F15.2,I6,F15.2,F9.3,F8.3)
      WRITE(6,707)SSW,IDW,SMW
707   FORMAT(/6X,'WITHIN',4X,F15.2,I6,F15.2)
      GO TO 711
710   WRITE(6,706)INAME,JNAME,KNAME,SSIJK,IDW,SMW
711   RETURN
      END
C**** WMU-AM: #1.9.3, MOD=1, MTO, 23-SEP-77 ****
C
C
C	PROPORTIONAL ANALYSIS
C	(WRITTEN BY DR. MICHAEL R. STOLINE, WMU, SEP-77)
C
C
C---------------ALL ARGS. ARE INPUT.
	SUBROUTINE PROP(IMAX,JMAX,KMAX,XI,XJ,XK,XIJ,XIK,XJK,
	1	XIJK,NI,NJ,NK,NIJ,NJK,NIK,NIJK,SSTOT,SSW,N,CT,SSC)
	DIMENSION XI(1),XJ(1),XK(1),XIJ(1),XIK(1),XJK(1),XIJK(1),
	1	NI(1),NJ(1),NK(1),NIJ(1),NJK(1),NIK(1),NIJK(1)
	SSI=0.
	SSJ=0.
	SSK=0.
	SSIJ=0.
	SSIK=0.
	SSJK=0.
	SSIJK=0.
	DO 1 I=1,IMAX
1	SSI=SSI+(XI(I)**2.)/NI(I)
	SSI=SSI-CT
	DO 2 J=1,JMAX
2	SSJ=SSJ+(XJ(J)**2.)/NJ(J)
	SSJ=SSJ-CT
	DO 3 K=1,KMAX
3	SSK=SSK+(XK(K)**2.)/NK(K)
	SSK=SSK-CT
	DO 4 I=1,IMAX
	DO 4 J=1,JMAX
	IJ=I+(J-1)*IMAX
4	SSIJ=SSIJ+(XIJ(IJ)**2.)/NIJ(IJ)
	DO 5 I=1,IMAX
	DO 5 K=1,KMAX
	IK=I+(K-1)*IMAX
5	SSIK=SSIK+(XIK(IK)**2.)/NIK(IK)
	DO 6 J=1,JMAX
	DO 6 K=1,KMAX
	JK=J+(K-1)*JMAX
6	SSJK=SSJK+(XJK(JK)**2.)/NJK(JK)
	SSIJ=SSIJ-CT
	SSIK=SSIK-CT
	SSJK=SSJK-CT
	SSIJ=SSIJ-SSI-SSJ
	SSIK=SSIK-SSI-SSK
	SSJK=SSJK-SSJ-SSK
	SSIJK=SSC-SSI-SSJ-SSK-SSIJ-SSIK-SSJK
	IDF=IMAX-1
	JDF=JMAX-1
	KDF=KMAX-1
	IJDF=IDF*JDF
	IKDF=IDF*KDF
	JKDF=JDF*KDF
	IJKDF=IJDF*KDF
	MSDF=N-IMAX*JMAX*KMAX
	NM1=N-1
	WRITE(6,101)
	WRITE(6,102)SSI,IDF,SSJ,JDF,SSK,KDF,SSIJ,IJDF,SSIK,IKDF,
	1	SSJK,JKDF,SSIJK,IJKDF,SSW,MSDF,SSTOT,NM1
101	FORMAT(/'0    THE SAMPLE SIZES ARE PROPORTIONAL.  THE SUMS'/
	1	' OF SQUARES ARE DECOMPOSABLE INTO MAIN EFFECTS SUMS'/
	1	' OF SQUARES AND INTERACTION SUMS OF SQUARES, AS IN'/
	1	' THE BALANCED CASE.  HOWEVER, THE USUAL F-TESTS ARE'/
	1	' NOT VALID IN THE PROPORTIONAL CASE, SO THESE TESTS'/
	1	' ARE NOT GIVEN IN THE PROPORTIONAL AOV TABLE.  A COM-'/
	1	' PLICATED EXACT ANALYSIS FOR PROPORTIONAL CASES IS'/
	1	' DESCRIBED IN BANCROFT, BUT IS NOT GIVEN HERE.'/
	1	'     A SUITABLE APPROXIMATE ANALYSIS FOR PROPORTIONAL'/
	1	' CASES IS THE LEAST SQUARES(FITTING CONSTANTS) METHOD'/
	1	' GIVEN HERE.'///
	1	20X,'PROPORTIONAL AOV TABLE'///
	1	' SOURCE',T15,'SUM OF SQUARES',T35,'DEGREES OF FREEDOM'/
	1	' ------',T15,'--------------',T35,'------------------')
102	FORMAT ('   A',T15,F11.3,T35,I11/
	1	'   B',T15,F11.3,T35,I11/
	1	'   C',T15,F11.3,T35,I11/
	1	'   AB',T15,F11.3,T35,I11/
	1	'   AC',T15,F11.3,T35,I11/
	1	'   BC',T15,F11.3,T35,I11/
	1	'   ABC',T15,F11.3,T35,I11/
	1	'   WITHIN',T15,F11.3,T35,I11/
	1	1X,51('-')/
	1	'   TOTAL',T15,F11.3,T35,I11//)

C**** END = PROP, NEW-ROUTINE ****

	END
