Trailing-Edge
-
PDP-10 Archives
-
decuslib10-10
-
43,50520/stp7.stp
There is 1 other file named stp7.stp in the archive. Click here to see a list.
C *** STAT PACK ***
C SUBROUTINE FOR 2 WAY ANALYSIS OF VARIANCE.
C CALLING SEQUENCE: CALL ANOV2(NV,NC,MV,MC,DATA,VMN,STD,NAMES,IX,IY)
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 IX - IS A VECTOR AT LEAST NC LONG
C IY - IS A VECTOR AT LEAST NC LONG
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,IX,IY)
DIMENSION DATA(MC,MV), VMN(1), STD(1), OPT(6), OPTD(6)
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),AM(20),YT(20,20),IBK(2),RNG1(3,20),RNG2(3,20)
DIMENSION NGRPS(2),C(20),IX(1),IY(1),LMOUT(15)
DIMENSION ICELLS(40),NAMES(1),ICLV(20,20),YR(20,20)
EQUIVALENCE (IBK(1),IBK1),(IBK(2),IBK2)
ISQ=5
IF(IOUT.EQ.21) ISQ=10
ALL=0
IF(ICC.NE.2) WRITE(IDLG,9910)
9910 FORMAT('0OPTION: GROUP IS NO LONGER AVAILABLE')
5 IF(ICC.NE.2) WRITE (IDLG,1)
OPTD(1)=0
OPTD(2)=0
OPTD(3)=0
OPTD(4)=0
OPTD(5)=0
OPTD(6)=0
1 FORMAT ('0LIST OPTIONS SEPARATED BY COMMAS')
READ (ICC,3) OPT
2 FORMAT ('0THE 2 WAY ANALYSIS OF VARIANCE ASSUMES CELLS TO '/
1' BE MADE UP OF INDIVIDUAL VARIABLES.'/
2' IT IS POSSIBLE HOWEVER TO BREAK 1 VARIABLE INTO'/
3' CELLS, BY MEANS OF TWO OTHER VARIABLES.'/
4' IN THIS CASE, RANGES FOR BREAK-DOWNS ARE ASSUMED'/
5' TO BE SUPPLIED BY THE USER. WHEN THE PROGRAM'/
6' CALLS FOR THE RANGES, THE USER MAY SPECIFY ANOTHER'/
7' OPTION WHICH ALLOWS THE PROGRAM TO AUTOMATICALLY FORM'/
8' THE BREAKDOWNS BY INDIVIDUAL VALUES OF THE')
632 FORMAT(' BREAKDOWN VARIABLE. IF MORE THAN 20'/
1' UNIQUE VALUES EXIST, THE PROGRAM AUTOMATICALLY'/
2' FORMS THE RANGES FOR BREAKDOWN'/ ' OPTIONS ARE:'/
3' "BARTL" - PERFORM BARTLETT''S TEST ON DATA.'/
5' "BREAK" - FORM GROUPINGS BASED ON VALUES OF BREAK VARIABLES.'/
6' "AUTO" - AUTOMATIC BREAKDOWN (SPECIFY WHEN ASKED FOR RANGES)'/
7' "DISCR" - BOTH GROUPS ARE AUTO (SPECIFY NOW)'/
8' "HEADR" - OMIT SIZE, MEANS, AND STD. DEV. REPORT'/
9' "RANGE" - DISPLAY BREAKDOWN RANGES IF AUTOMATICALLY FORMED'/
2'0IF NO OPTIONS ARE DESIRED TYPE A RETURN')
3 FORMAT (6(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 90
OPTD(1)=1
GO TO 4
90 IF(OPT(I).NE.'WMEAN') GO TO 7
WRITE(IDLG,9001)
9001 FORMAT(' "WMEAN" IS NOW AUTOMATIC - DO NOT USE THIS OPTION')
GO TO 4
7 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 550
OPTD(5)=1
GO TO 4
550 IF(OPT(I).NE.'BARTL') GO TO 509
OPTD(6)=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
NFAC1='ONE'
NFAC2='TWO'
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 513
WRITE (IDLG,14)
14 FORMAT ('0MAXIMUM OF 20')
GO TO 11
513 DO 514 I=1,NGRPA
514 ENCODE(5,515,RNG1(3,I)) I
515 FORMAT(I4,1X)
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 516
WRITE (IDLG,14)
GO TO 15
516 DO 517 I=1,NGRPB
517 ENCODE(5,515,RNG2(3,I)) I
18 IF(ICC.NE.2) WRITE (IDLG,19)
19 FORMAT ('0TYPE IN EACH VARIABLE AFTER THE'/
1' CORRESPONDING LEVEL, "EMPTY"-INDICATES AN EMPTY CELL')
XT=NC
L=0
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(ICLV(I,J),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(ICLV(I,J).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 ICLV(I,J)=L
GO TO 26
441 IF(ICLV(I,J).EQ.0)GO TO 26
IF(I.EQ.1) GO TO 444
DO 442 K=1,I-1
DO 442 II=1,NGRPB
IF(ICLV(K,II).EQ.ICLV(I,J)) GO TO 24
442 CONTINUE
444 IF(J.EQ.1) GO TO 26
DO 440 K=1,J-1
IF(ICLV(I,K).EQ.ICLV(I,J)) GO TO 24
440 CONTINUE
GO TO 26
24 WRITE(IDLG,443)
443 FORMAT(' THIS VARIABLE USED PREVIOUSLY')
GO TO 22
26 CONTINUE
IF(L.EQ.0) GO TO 407
L=L-1
DO 410 I=1,NGRPA
DO 410 J=1,NGRPB
IF(ICLV(I,J).GE.0) GO TO 410
ICLV(I,J)=L-ICLV(I,J)
410 CONTINUE
400 M=0
DO 402 I=1,NGRPA
DO 402 J=1,NGRPB
IF(ICLV(I,J).GE.0) GO TO 402
ICLV(I,J)=ICLV(I,J)-1
IF(ICLV(I,J).GE.-(NV-M)) GO TO 403
M=M+1
402 CONTINUE
RETURN
403 IF(M.EQ.0) GO TO 405
N=0
DO 404 K=1,NGRPA
DO 404 L=1,NGRPB
IF(ICLV(K,L).GE.0) GO TO 404
ICLV(K,L)=ICLV(I,J)-M+N
N=N+1
IF(N.EQ.M) GO TO 405
404 CONTINUE
405 DO 406 I=1,NGRPA
DO 406 J=1,NGRPB
DO 406 K=1,NGRPA
DO 406 L=1,NGRPB
IF((I.EQ.K).AND.(J.EQ.L)) GO TO 406
IF(IABS(ICLV(I,J)).EQ.IABS(ICLV(K,L))) GO TO 400
406 CONTINUE
407 DO 88 I=1,20
DO 88 J=1,20
T(I,J)=0
YT(I,J)=0
88 YS(I,J)=0
DO 25 I=1,NGRPA
DO 25 J=1,NGRPB
ICELL=ICLV(I,J)
IF(ICELL.EQ.0) GO TO 25
IF(ICELL.LT.0) ICELL=-ICELL
C YT(I,J)=VMN(ICELL)*XT
C T(I,J)=XT
C YS(I,J)=((STD(ICELL)**2)*(XT*(XT-1.))+YT(I,J)**2)/XT
C** IT COSTS MORE, BUT WHAT THE HECK, IT GETS BETTER ANSWERS
T(I,J)=XT
DO 2611 K=1,NC
YT(I,J)=YT(I,J)+DATA(K,ICELL)
2611 YS(I,J)=YS(I,J)+DATA(K,ICELL)**2
25 CONTINUE
GO TO 87
C
C CELLS ON INDIV. VAR. BY BREAK-DOWN.
C
30 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 30
IF(IERR.EQ.1) GO TO 30
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(IBK(I),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(IBK(I).GT.0) GO TO 40
WRITE(IDLG,450)
450 FORMAT(' *, ?, ALL MAY NOT BE USED HERE')
GO TO 42
40 IDIS1=OPTD(4)
NFAC1=NAMES(IBK1)
CALL RANGES(NV,NC,MV,MC,DATA,IDIS1,RNG1,20,NGRPA,IRET,
1IBK1,NAMES,IX)
IF(IRET.EQ.1) RETURN
IDIS2=OPTD(4)
NFAC2=NAMES(IBK2)
CALL RANGES(NV,NC,MV,MC,DATA,IDIS2,RNG2,20,NGRPB,IRET,IBK2,
1NAMES,IX)
IF(IRET.EQ.1) RETURN
IF((OPTD(5).NE.1).OR.(IDIS1.NE.1)) GO TO 50
WRITE(IDLG,47) NAMES(IBK1)
47 FORMAT('0BREAKDOWN RANGES FOR VARIABLE: ',A5)
DO 48 I=1,NGRPA
48 WRITE(IDLG,49)(RNG1(J,I),J=1,2)
49 FORMAT(1X,G10.4,' , ',G10.4)
50 IF((OPTD(5).NE.1).OR.(IDIS2.NE.1)) GO TO 52
WRITE(IDLG,47) NAMES(IBK2)
DO 51 I=1,NGRPA
51 WRITE(IDLG,49)(RNG2(J,I),J=1,2)
52 DO 43 I=1,NC
DO 44 J=1,NGRPA
IF(DATA(I,IBK1).LT.RNG1(1,J)) GO TO 44
IF(DATA(I,IBK1).GT.RNG1(2,J)) GO TO 44
IX(I)=J
GO TO 45
44 CONTINUE
IX(I)=0
GO TO 43
45 DO 46 J=1,NGRPB
IF(DATA(I,IBK2).LT.RNG2(1,J)) GO TO 46
IF(DATA(I,IBK2).GT.RNG2(2,J)) GO TO 46
IY(I)=J
GO TO 43
46 CONTINUE
IY(I)=0
43 CONTINUE
C
C AT THIS POINT RANGES HAVE BEEN DET.
C BREAK DOWN VARIABLES
C
KCELL=0
2000 KCELL=KCELL+1
IF(KCELL.GT.NZZ) RETURN
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,NGRPA
DO 388 J=1,NGRPB
T(I,J)=0
YT(I,J)=0
388 YS(I,J)=0
DO 60 I=1,NC
J=IX(I)
K=IY(I)
IF(J.EQ.0) GO TO 60
IF(K.EQ.0) GO TO 60
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.
60 CONTINUE
C
C BEGINNING OF THE ANALYSIS OF VARIANCE PORTION
C T(I,J) SHOULD CONTAIN THE FREQUENCY OF THE CELL I,J
C YS(I,J) SHOULD CONTAIN THE SUM OF X'S SQUARED FOR CELL I,J
C YT(I,J) SHOULD CONTAIN THE SUM OF THE X'S FOR CELL I,J
C
87 DO 91 I=1,20
S1(I)=0
S2(I)=0
91 CONTINUE
C
C FORM THE N'S ON LEVELS FOR BOTH FACTORS
C
DO 65 I=1,NGRPA
DO 66 J=1,NGRPB
IF(T(I,J).EQ.0) GO TO 66
S1(I)=S1(I)+T(I,J)
S2(J)=S2(J)+T(I,J)
66 CONTINUE
65 CONTINUE
C
C NOW CHECK TO SEE THAT ALL LEVELS HAVE AT LEAST 1 OBSERVATION
C
70 N1=NGRPA
N2=NGRPB
IF(N1.GT.1) GO TO 127
WRITE(IDLG,126) NFAC1
126 FORMAT(' FEWER THAN 2 LEVELS FOR FACTOR ',A5)
GO TO 1000
127 IF(N2.GT.1) GO TO 131
WRITE(IDLG,126) NFAC2
GO TO 1000
131 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 AND WILL BE ELIMINATED')
75 WRITE(IDLG,76) NFAC1,RNG1(3,I)
76 FORMAT(' FACTOR ',A5,' LEVEL: ',A5)
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,76) NFAC2,RNG2(3,J)
77 CONTINUE
IF(SWTCH.EQ.0) GO TO 80
C
C FOLLOWING USED ONLY IF LEVELS ARE ELIMINATED
C
L=1
DO 82 I=1,NGRPA
IF(S1(I).EQ.0) GO TO 82
IF(I.EQ.L) GO TO 81
DO 83 J=1,NGRPB
YT(L,J)=YT(I,J)
YS(L,J)=YS(I,J)
T(L,J)=T(I,J)
83 CONTINUE
81 L=L+1
82 CONTINUE
NGRPA=L-1
L=1
DO 85 J=1,NGRPB
IF(S2(J).EQ.0) GO TO 85
IF(J.EQ.L) GO TO 84
DO 86 I=1,NGRPA
YT(I,L)=YT(I,J)
YS(I,L)=YS(I,J)
T(I,L)=T(I,J)
86 CONTINUE
84 L=L+1
85 CONTINUE
NGRPB=L-1
GO TO 87
C
C AT THIS POINT EACH LEVEL HAS AT LEAST ONE OBSERVATION
C SO EVERYTHING OK TO CONTINUE
C
80 DO 540 I=1,N1
Y1(I)=0
540 AM(I)=0
DO 541 J=1,N2
C(J)=0
B1(J)=0
541 Y2(J)=0
RAB=0
SST=0
SS=0
G=0
SSC=0
SSAB=0
AN=0
SSBA=0
TW=0
TC=0
IWMEAN=0
C
C FORM REMAINING TOTALS FOR THE ANALYSIS
C
DO 542 I=1,NGRPA
DO 543 J=1,NGRPB
IF(T(I,J).EQ.0) GO TO 544
TC=TC+1.
TW=TW+T(I,J)-1
Y1(I)=Y1(I)+YT(I,J)
Y2(J)=Y2(J)+YT(I,J)
SSC=SSC+YT(I,J)**2/T(I,J)
SST=SST+YS(I,J)
YS(I,J)=YS(I,J)-(YT(I,J)**2)/T(I,J)
IF(T(I,J).EQ.1) GO TO 543
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)
GO TO 543
544 IWMEAN=1
543 CONTINUE
542 CONTINUE
DO 545 J=1,NGRPB
545 SSBA=SSBA+(Y2(J)**2)/S2(J)
DO 93 I=1,N1
SSAB=SSAB+Y1(I)**2/S1(I)
SS=SS+S1(I)
93 G=G+Y1(I)
G1=(G*G)/SS
SST=SST-G1
SSC=SSC-G1
SSAB=SSAB-G1
SSBA=SSBA-G1
SSW=SST-SSC
SMC=SSC/(TC-1.)
IF(TW.EQ.0) GO TO 505
SMW=SSW/TW
IF(SMW.EQ.0) GO TO 505
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 551
C
C MEANS AND STANDARD DEVIATION REPORT
C
DO 116 I=1,N2,ISQ
IF(IOUT.NE.21) GO TO 324
LINES=LINES+4
IF(LINES.LE.(LINPP-8)) GO TO 324
CALL PRNTHD
LINES=6
324 WRITE (IOUT,107) NFAC1,NFAC2
107 FORMAT('0FACTOR'/1X,A5,25X,'FACTOR ',A5)
LINES=LINES+4
NEND=I+ISQ-1
IF(NEND.GT.N2) NEND=N2
WRITE(IOUT,108)(RNG2(3,J),J=I,NEND)
108 FORMAT(' LEVEL',5X,10(3X,A5,4X))
DO 322 K=1,N1
IF(IOUT.NE.21) GO TO 325
LINES=LINES+4
IF(OPTD(1).NE.1) LINES=LINES+1
IF(LINES.LE.LINPP) GO TO 325
CALL PRNTHD
WRITE(IOUT,107) NFAC1,NFAC2
WRITE(IOUT,108)(RNG2(3,J),J=I,NEND)
LINES=10
IF(OPTD(1).NE.1) LINES=11
325 DO 326 J=I,NEND
326 LMOUT(J-I+1)=T(K,J)
LM=NEND-I+1
IF(OPTD(1).EQ.1) GO TO 115
WRITE(IOUT,519) RNG1(3,K),(NAMES(IABS(ICLV(K,J))),J=I,NEND)
519 FORMAT('0',A5,5X,10(3X,A5,4X))
WRITE(IOUT,520)(LMOUT(J),J=1,LM)
520 FORMAT(9X,'N',1X,10(I8,4X))
GO TO 521
115 WRITE(IOUT,110) RNG1(3,K),(LMOUT(J),J=1,LM)
110 FORMAT('0',A5,3X,'N',1X,12(I8,4X))
521 WRITE(IOUT,114)(YT(K,J),J=I,NEND)
114 FORMAT(6X,'MEAN',2X,10(G11.5,1X))
WRITE(IOUT,154)(YS(K,J),J=I,NEND)
154 FORMAT(5X,'STDEV',2X,10(G11.5,1X))
322 CONTINUE
IF(IOUT.NE.21) GO TO 111
LINES=LINES+2
IF(LINES.GT.LINPP) GO TO 116
111 WRITE(IOUT,113)
113 FORMAT(//)
116 CONTINUE
C
C BARTLETT'S TEST
C
551 IF(OPTD(6).EQ.0) GO TO 599
A=0
SUMFS=0
SUM1F=0
SUMF=0
SUMLS=0
DO 552 I=1,NGRPA
DO 552 J=1,NGRPB
DF=T(I,J)-1.
IF(DF.LE.0) GO TO 552
SUMF=SUMF+DF
SUM1F=SUM1F+1./DF
SUMLS=SUMLS+ALOG(YS(I,J)**2)*DF
SUMFS=SUMFS+DF*YS(I,J)**2
A=A+1
552 CONTINUE
SBAR=SUMFS/SUMF
AMM=SUMF*ALOG(SBAR)-SUMLS
CC=1.+1./(3.*(A-1.))*(SUM1F-(1./SUMF))
CHI=AMM/CC
KK=A-1.
IF(KK.LE.0) GO TO 599
IF(IOUT.NE.21) GO TO 553
LINES=LINES+5
IF(LINES.LE.LINPP) GO TO 553
CALL PRNTHD
LINES=7
553 WRITE(IOUT,554) CHI,KK
554 FORMAT('0BARTLETT''S TEST OF HOMOGENEITY OF VARIANCE'/
1' CHI SQUARE =',G12.4/' WITH ',I3,' DEGREES OF FREEDOM')
C**RRB/13OCT77
IERR=0
CALL CHIPRB(CHI,KK,PROB)
IF(PROB.EQ.99)IERR=1
IF(IERR.EQ.0) GO TO 556
WRITE(IOUT,557)
557 FORMAT(' NO PROBABILITY CALCULATED DF OR CHI OUTSIDE LIMITS')
GO TO 599
556 WRITE(IOUT,558) PROB
558 FORMAT(' HAVING A PROBABILITY OF ',F6.3)
C
C PRELIMINARY ANOVA (ELIMINATED IF TW=0 OR IF BALANCED ANOVA)
C
599 IF(TW.EQ.0) GO TO 335
DO 330 I=1,NGRPA
DO 330 K=1,NGRPB
IF(T(I,K).NE.T(1,1)) GO TO 155
330 CONTINUE
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(9X,'SOURCE',14X,'DF',8X,'SS',8X,'MS',8X,'F',4X,'PROB')
KK=TC-1.
NDG=TW
IF(SMW.NE.0) PROB=FISHER(KK,NDG,F)
IF(SMW.NE.0) WRITE (IOUT,119)KK,SSC,SMC,F,PROB
IF(SMW.EQ.0) WRITE(IOUT,119) KK,SSC,SMC
119 FORMAT(9X'CELLS',12X,I5,3F10.2,2X,F6.4)
KK=N1-1
WRITE (IOUT,120)NFAC1,NFAC2,KK,SSAB
120 FORMAT(2X,A5,' IGNORING ',A5,4X,I5,F10.2)
KK=N2-1
WRITE (IOUT,120)NFAC2,NFAC1,KK,SSBA
KK=TW
WRITE (IOUT,121)KK,SSW,SMW
121 FORMAT(9X,'WITHIN',11X,I5,2F10.2)
KK=SS-1.
WRITE (IOUT,122)KK,SST
122 FORMAT(9X,'TOTAL',12X,I5,F10.2)
C
C DO CALCULATIONS FOR THE LEAST SQUARES ANOVA
C
335 DO 125 I=1,N2
DO 125 J=1,N2
YR(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
128 YR(I,J)=YR(I,J)+(T(K,I)*T(K,J))/S1(K)
DO 129 I=1,N2
129 YR(I,I)=YR(I,I)-S2(I)
DO 133 J=1,N2
130 DO 132 K=1,N1
132 C(J)=C(J)+(T(K,J)*Y1(K))/S1(K)
YR(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(((YR(I,I)+100.)-100.).NE.0.0) GO TO 220
IF(I.EQ.N2) GO TO 260
DO 211 J=I+1,N2
IF(((YR(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
YR(I,K)=YR(I,K)+YR(J,K)
213 YS(I,K)=YS(I,K)+YS(J,K)
220 GOB=YR(I,I)
DO 221 J=1,N2
YR(I,J)=YR(I,J)/GOB
221 YS(I,J)=YS(I,J)/GOB
DO 230 L=1,N2
IF(L.EQ.I) GO TO 230
GOB=YR(L,I)
DO 231 J=1,N2
YR(L,J)=YR(L,J)-GOB*YR(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 CONTINUE
C 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.
IF(TI.NE.0) SMI=SI/TI
TAB=RAB-SSBA
TMAB=TAB/(P1-1.)
TBA=RAB-SSAB
TMBA=TBA/(P2-1.)
NDG=TW
IF(TW.NE.0) GO TO 105
NDG=TI
SMW=SMI
105 IF(SMW.EQ.0) GO TO 144
FI=SMI/SMW
FA=TMAB/SMW
FB=TMBA/SMW
C
C LEAST SQUARES ANOVA OUTPUT
C
144 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(//29X,'LEAST SQUARES ANOVA')
WRITE (IOUT,118)
KK=TC-1.
WRITE (IOUT,119)KK,SSC
KK=P1-1.
IF(SMW.NE.0) PROB=FISHER(KK,NDG,FA)
IF(SMW.NE.0) WRITE (IOUT,146) NFAC1,NFAC2,KK,TAB,TMAB,FA,PROB
IF(SMW.EQ.0) WRITE(IOUT,146) NFAC1,NFAC2,KK,TAB,TMAB
146 FORMAT(1X,A5,' ELIMINATING ',A5,2X,I5,3F10.2,2X,F6.4)
KK=P2-1.
IF(SMW.NE.0) PROB=FISHER(KK,NDG,FB)
IF(SMW.NE.0) WRITE (IOUT,146) NFAC2,NFAC1,KK,TBA,TMBA,FB,PROB
IF(SMW.EQ.0) WRITE(IOUT,146) NFAC2,NFAC1,KK,TBA,TMBA
KK=TI
IF(TW.NE.0) PROB=FISHER(KK,NDG,FI)
IF(TW.NE.0) WRITE (IOUT,148) NFAC1,NFAC2,KK,SI,SMI,FI,PROB
IF(TW.EQ.0) WRITE(IOUT,148) NFAC1,NFAC2,KK,SI,SMI
148 FORMAT(5X,A5,' BY ',A5,7X,I5,3F10.2,2X,F6.4)
KK=TW
IF(TW.NE.0) WRITE (IOUT,121)KK,SSW,SMW
KK=SS-1
WRITE(IOUT,122) KK,SST
C
C WEIGHTED MEANS ANALYSIS OF VARIANCE
C
IF(IWMEAN.EQ.0)GO TO 808
WRITE(IOUT,807)
807 FORMAT(//,' WEIGHTED MEANS ANALYSIS CANNOT BE PERFORMED',
1' WITH EMPTY CELLS')
GO TO 1000
C
808 IF(IOUT.NE.21)GO TO 1842
IF(LINES+9.LE.LINPP)GO TO 1841
CALL PRNTHD
LINES=2
1841 LINES=LINES+9
1842 WRITE(IOUT,1840)
1840 FORMAT(/,
1' THE FOLLOWING OUTPUT IS USED FOR ANALYSIS OF UNBALANCED',/,
2' TWO-WAY ANALYSIS OF VARIANCE DESIGNS, WHICH IS OUTLINED',/,
3' IN SECTION 1.9 OF T.A. BANCROFT''S "TOPICS IN INTERMED-',/,
4' IATE STATISTICAL METHODS".',/,
5' CONSULT THE DOCUMENTATION OF THE WMU COMPUTER CENTER',/,
6' PROGRAM #1.9.2(TWOAOV) FOR A DETAILED DESCRIPTION OF THE',/,
7' CORRECT USE OF THIS OUTPUT.',/)
C
C CHECK FOR BALANCED CASE
C
DO 809 I=1,N1
DO 809 J=1,N2
IF(T(I,J).NE.T(1,1)) GO TO 810
809 CONTINUE
WRITE(IOUT,1822)
1822 FORMAT(/,' THE DESIGN IS BALANCED. IN THIS CASE THE WEIGHTED',/,
#' MEANS AND LEAST SQUARES ANOVA''S ARE IDENTICAL.',//)
GO TO 1000
C
C CHECK FOR PROPORTIONALITY
C
1823 DO 1824 I=2,N1
TI1=T(I,1)
DO 1824 J=2,N2
IF(TI1*T(I,J).NE.T(1,I)*TI1)GO TO 810
1824 CONTINUE
IF(IOUT.NE.21)GO TO 1825
IF(LINES+9.LE.LINPP)GO TO 1826
CALL PRNTHD
1826 LINES=0
1825 LINES=LINES+9
WRITE(IOUT,1831)
1831 FORMAT(/,' THE DESIGN IS PROPORTIONAL, BUT NOT BALANCED.',/,
#' IN THIS CASE THE LEAST SQUARES ANOVA GIVES THE',/,
#' CORRECT SUM OF SQUARES AND MEAN SQUARES, BUT NOT',/,
#' CORRECT F-TESTS. SEE T.A. BANCROFT "TOPICS IN',/,
#' INTERMEDIATE STATISTICAL METHODS(P. 14-15)" FOR',/,
#' THE PROPER F-TEST PROCEDURES FOR THE PROPOR-',/,
#'TIONAL CASE.',//)
GO TO 1000
810 DO 811 I=1,N1
Y1(I)=0
811 S1(I)=0
DO 812 J=1,N2
Y2(J)=0
812 S2(J)=0
DO 813 I=1,N1
DO 814 J=1,N2
Y1(I)=Y1(I)+YT(I,J)
S1(I)=S1(I)+1.0/T(I,J)
Y2(J)=Y2(J)+YT(I,J)
814 S2(J)=S2(J)+1.0/T(I,J)
Y1(I)=Y1(I)/N1
813 S1(I)=N1**2/S1(I)
DO 815 J=1,N2
Y2(J)=Y2(J)/N2
815 S2(J)=N2**2/S2(J)
A=0
B=0
D=0
DO 816 I=1,N1
A=A+S1(I)*Y1(I)**2
B=B+S1(I)*Y1(I)
816 D=D+S1(I)
SSR=A-B**2/D
A=0
B=0
D=0
DO 817 J=1,N2
A=A+S2(J)*Y2(J)**2
B=B+S2(J)*Y2(J)
817 D=D+S2(J)
SSCO=A-B**2/D
SMR=SSR/(N1-1.)
SMCO=SSCO/(N2-1.)
C
C OUTPUT OF WEIGHTED MEANS ANOVA
C
IF(IOUT.NE.21) GO TO 818
LINES=LINES+10
IF(LINES.LE.LINPP) GO TO 818
CALL PRNTHD
LINES=10
818 WRITE(IOUT,822)
822 FORMAT(//20X,'WEIGHTED MEANS ANOVA')
WRITE(IOUT,118)
KK=TC-1
WRITE(IOUT,119) KK,SSC,SMC
KK=P1-1
F=SMR/SMW
PROB=FISHER(KK,NDG,F)
WRITE(IOUT,819) NFAC1,KK,SSR,SMR,F,PROB
819 FORMAT(10X,A5,11X,I5,3F10.2,2X,F6.4)
KK=P2-1
F=SMCO/SMW
PROB=FISHER(KK,NDG,F)
WRITE(IOUT,819) NFAC2,KK,SSCO,SMCO,F,PROB
KK=TI
IF(TW.EQ.0) GO TO 820
PROB=FISHER(KK,NDG,FI)
WRITE(IOUT,148) NFAC1,NFAC2,KK,SI,SMI,FI,PROB
KK=TW
WRITE(IOUT,121) KK,SSW,SMW
GO TO 821
820 WRITE(IOUT,148) NFAC1,NFAC2,TI,SI,SMI
821 KK=SS-1
WRITE(IOUT,122) KK,SST
1000 IF(OPTD(1).NE.1) GO TO 400
GO TO 2000
END
SUBROUTINE RANGES(NV,NC,MV,MC,DATA,IDISCR,RANGE,MAXRN,NRANG,
1IRET,IVAR,NAMES,D)
DIMENSION DATA(MC,MV),RANGE(3,1),NAMES(1),IN(40),LDAT(15),MDAT(3)
DIMENSION D(1)
COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
DO 40 I=1,NC
40 D(I)=DATA(I,IVAR)
IRET=0
IF(IDISCR.EQ.1) GO TO 50
IF(ICC.NE.2) WRITE(IDLG,1) NAMES(IVAR)
1 FORMAT(' ENTER RANGES FOR VARIABLE ',A5/)
NRANG=1
2 IF(ICC.NE.2) WRITE(IDLG,6) NRANG
6 FORMAT('+',I2,'? ',$)
READ(ICC,3,END=32) IN
3 FORMAT(40A1)
IF((IN(1).EQ.' ').AND.(IN(2).EQ.' ').AND.(IN(3).EQ.' ')) GO TO 32
IF((IN(1).EQ.'S').AND.(IN(2).EQ.'T').AND.(IN(3).EQ.'O').AND.
1(IN(4).EQ.'P')) GO TO 32
IF((IN(1).EQ.'A').AND.(IN(2).EQ.'U').AND.(IN(3).EQ.'T').AND.
1(IN(4).EQ.'O')) GO TO 50
IF(IN(1).NE.'!') GO TO 4
IRET=1
RETURN
4 IF((IN(1).NE.'H').OR.(IN(2).NE.'E').OR.(IN(3).NE.'L').OR.
1(IN(4).NE.'P')) GO TO 7
WRITE(IDLG,5) MAXRN,MAXRN
5 FORMAT('0RANGES ARE ENTERED ONE PER LINE IN RESPONSE TO A'/
1' QUESTION MARK. UP TO ',I2,' RANGES MAY BE ENTERED, EACH RANGE'/
2' CONSISTING OF THE MINIMUM FOR THE RANGE, A COMMA, THE'/
3' MAXIMUM FOR THE RANGE, AND IF DESIRED A COMMA AND THE NAME'/
4' FOR THE RANGE. WHEN THE LAST RANGE HAS BEEN ENTERED, TYPE A'/
5' BLANK LINE, "STOP", OR CONTROL Z(^Z).'/
6'0RANGES MAY BE AUTOMATICALLY CREATED BY SPECIFYING "AUTO/", AND'/
7' THE NUMBER OF RANGES THE DATA IS TO BE BROKEN INTO. IF THE'/
8' "/" AND NUMBER OF RANGES ARE NOT ENTERED, THE MAXIMUM (',I2,')'/
9' NUMBER OF RANGES WILL BE ASSUMED. IF THERE EXISTS FEWER'/
1' INDIVIDUAL VALUES THAN THE NUMBER OF RANGES SPECIFIED, A RANGE'/
2' WILL BE CREATED FOR EACH INDIVIDUAL VALUE. IF MORE VALUES'/
3' EXIST THAN RANGES, THEN THE RANGES WILL BE CREATED BY BREAKING'/
4' THE DIFFERENCE BETWEEN THE MAXIMUM AND MINIMUM VALUES INTO'/
5' EQUAL SIZED RANGES. RANGES NOT CONTAINING ANY VALUES WILL BE'/
6' ELIMINATED.'/)
GO TO 2
7 M=1
I=0
23 L=1
IDP=0
IEXP=0
MINUS=0
DO 8 J=1,15
8 LDAT(J)=' '
I=I+1
9 IF(IN(I).EQ.',') GO TO 20
IF(IN(I).EQ.' ') GO TO 20
IF(IN(I).EQ.'-') GO TO 18
IF(IN(I).EQ.'E') GO TO 16
IF(IN(I).EQ.'.') GO TO 14
IF((IN(I).GE.'0').AND.(IN(I).LE.'9')) GO TO 11
WRITE(IDLG,10) IN(I)
10 FORMAT(' ILLEGAL CHARACTER (',A1,') IN ONE VALUE OF THE RANGE'/)
GO TO 2
11 IF(L.LE.15) GO TO 13
WRITE(IDLG,12)
12 FORMAT(' NUMBER GREATER THAN 15 CHARACTERS'/)
GO TO 2
13 LDAT(L)=IN(I)
L=L+1
I=I+1
GO TO 9
14 IDP=IDP+1
IF(IDP.LT.2) GO TO 13
WRITE(IDLG,15)
15 FORMAT(' ONLY 1 DECIMLE POINT PER NUMBER'/)
GO TO 2
16 IEXP=IEXP+1
MINUS=0
IF(IEXP.LT.2) GO TO 13
WRITE(IDLG,17)
17 FORMAT(' ONLY 1 EXPONENT PER NUMBER'/)
GO TO 2
18 MINUS=MINUS+1
IF(MINUS.LT.2) GO TO 13
WRITE(IDLG,19)
19 FORMAT(' ONLY 1 MINUS PER NUMBER'/)
GO TO 2
20 IF(LDAT(15).NE.' ') GO TO 39
DO 38 J=14,1,-1
38 LDAT(J+1)=LDAT(J)
LDAT(1)=' '
GO TO 20
39 ENCODE(15,3,MDAT) LDAT
DECODE(15,21,MDAT) RANGE(M,NRANG)
21 FORMAT(G15.0)
M=M+1
IF(M.EQ.3) GO TO 24
IF(IN(I).EQ.',') GO TO 23
WRITE(IDLG,22)
22 FORMAT(' NO COMMA SEPARATING MINIMUM FORM MAXIMUM IN RANGE'/)
GO TO 2
24 IF(RANGE(1,NRANG).LE.RANGE(2,NRANG)) GO TO 34
WRITE(IDLG,33)
33 FORMAT(' RANGES MUST BE ENTERED MINIMUM FIRST'/)
GO TO 2
34 IF(IN(I).EQ.' ') GO TO 28
C
C NAME
C
DO 25 J=1,5
25 LDAT(J)=' '
L=1
26 I=I+1
IF(IN(I).EQ.' ') GO TO 27
LDAT(L)=IN(I)
L=L+1
IF(L.LE.5) GO TO 26
27 IF(L.EQ.1) GO TO 28
ENCODE(5,3,RANGE(3,NRANG))(LDAT(L),L=1,5)
GO TO 30
28 ENCODE(5,29,RANGE(3,NRANG)) NRANG
29 FORMAT(2X,I2,1X)
30 IF(NRANG.LT.2) GO TO 37
DO 35 I=1,NRANG-1
IF(RANGE(1,I).GT.RANGE(2,NRANG)) GO TO 35
IF(RANGE(2,I).LT.RANGE(1,NRANG)) GO TO 35
WRITE(IDLG,36) RANGE(3,I)
36 FORMAT(' THE RANGE OVERLAPS RANGE ',A5/)
GO TO 2
35 CONTINUE
37 DO 41 I=1,NC
IF(D(I).LT.RANGE(1,NRANG))GO TO 41
IF(D(I).GT.RANGE(2,NRANG)) GO TO 41
GO TO 43
41 CONTINUE
WRITE(IDLG,42)
42 FORMAT(' NO OBSERVATIONS IN RANGE'/)
GO TO 2
43 NRANG=NRANG+1
IF(NRANG.LE.MAXRN) GO TO 2
WRITE(IDLG,31) MAXRN
31 FORMAT(' NO MORE RANGES ACCEPTED MAXIMUM (',I2,
1') ALREADY ENTERED'/)
32 NRANG=NRANG-1
RETURN
C
C AUTO WAS USED
C
50 NRANG=MAXRN
IF(IN(5).NE.'/') GO TO 60
IF(IDISCR.EQ.1) GO TO 60
I=6
51 IF((IN(I).GT.'9').OR.(IN(I).LT.'0')) GO TO 54
I=I+1
IF(I.LE.7) GO TO 51
IF(IN(I).EQ.' ') GO TO 54
52 WRITE(IDLG,53) MAXRN
53 FORMAT(' MAXIMUM NUMBER OF RANGES POSSIBLE IS ',I2/)
GO TO 2
54 IF(IN(7).NE.' ') GO TO 55
IN(7)=IN(6)
IN(6)=' '
55 ENCODE(2,3,MDAT)(IN(J),J=6,7)
DECODE(2,56,MDAT) NRANG
56 FORMAT(I2)
IF(NRANG.GT.MAXRN) GO TO 52
C
C BEGIN CREATION OF RANGES USE BINARY TABLE LOOKUP TO CONSTRUCT
C FREQUENCTY TABLE, AT SAME TIME ALSO CHECK FOR MIN AND MAX
C IF TO0 MANY INDIVIDUAL VALUES EXIST THE THESE WILL BE USED
C
60 AMAX=D(1)
AMIN=AMAX
RANGE(1,1)=AMAX
L=1
DO 61 I=2,NC
X=D(I)
IF(X.GT.AMAX) AMAX=X
IF(X.LT.AMIN) AMIN=X
IF(L.GT.NRANG) GO TO 61
IM=0
IF(X.LT.RANGE(1,1)) GO TO 64
IF(X.EQ.RANGE(1,1)) GO TO 61
IF(X.GT.RANGE(1,L)) GO TO 66
IF(X.EQ.RANGE(1,L)) GO TO 61
I1=1
I2=L
62 IM=(I1+I2)/2
IF(IM.EQ.I1) GO TO 64
IF(X.EQ.RANGE(1,IM)) GO TO 61
IF(X.LT.RANGE(1,IM)) GO TO 63
I1=IM
GO TO 62
63 I2=IM
GO TO 62
64 L=L+1
IF(L.GT.NRANG) GO TO 61
DO 65 J=L,IM+2,-1
65 RANGE(1,J)=RANGE(1,J-1)
RANGE(1,IM+1)=X
GO TO 61
66 L=L+1
IF(L.GT.NRANG) GO TO 61
RANGE(1,L)=X
61 CONTINUE
IF(L.GT.NRANG) GO TO 70
NRANG=L
DO 68 I=1,NRANG
RANGE(2,I)=RANGE(1,I)
ENCODE(5,29,RANGE(3,I)) I
68 CONTINUE
GO TO 72
70 RNG=(AMAX-AMIN)/NRANG
RANGE(1,1)=AMIN
DO 71 I=1,NRANG-1
RANGE(2,I)=RANGE(1,I)+RNG
RANGE(1,I+1)=RANGE(2,I)
ENCODE(5,29,RANGE(3,I)) I
71 CONTINUE
RANGE(2,NRANG)=AMAX
ENCODE(5,29,RANGE(3,NRANG)) NRANG
72 IDISCR=1
RETURN
END