C THIS SUBROUTINE WAS DESIGNED BY PROF. J. OGILVIE UNIVERSITY OF TORONTO C C YOU HAVE TO WRITE A MAINLINE PROGRAM TO GET THE DATA INTO ANOVAP C AND IT WILL DO THE REST C It is substantually more accurate than many other similar routines C under some circumstances as it avoids take the difference of very C large number to obtain a small difference. C C PURPOSE: C This subroutine will perform a factorial analysis of variance C for a 1 way to 7 way classification with one or more C observations per cell. C Means for main effects, 2 way and 3 way interactions as specified C and the analysis of variance table are printed. C The observations are replaced by residuals. C THIS IS VERY HANDY FOR MISSING DATA REPLACEMENT BY THE METHOD C OF RESIDUAL REDUCTION. JUST REPLACE THE MISSING VALUSE WITH C ANYTHING (A MEAN OR ZERO WILL DO) THEN RUN ANOVAP. TAKE C THE RESIDUAL THAT IS RETURNED IN PLACE OF THE VALUE YOU C SUPPLIED AND SUBTRACT IT FROM THE FIRST GUESS (THE MEAN C OR ZERO MENTIONED ABOVE, WATCH YOUR SIGNS) THEN USE THIS C AS A NEW ESTIMATE AND RUN THE PROGRAM AGAIN. CONTINUE C UNTILL THE RESIDUAL IS ARBITRARILY SMALL. REDUCE THE ERROR C DEGREES OF FREEDOM THE NUMBER OF MISSING DATA POINTS. C (See "Missing Values in Experiments Analysed on Automatic C Computers" Applied Statistics vol. 5, pp203-6 1956) C C HOW TO USE: C CALLANOVAP(X,N4,NCLASS,NCELL,LEVEL,VARNAM,MEANS,T,TTT,MDF,NPOOL) C C where: C C X C is a floating point vector with the data arranged in a specific order, C described below. C After ANOVAP is completed, the X is replaced by residuals. C X must be dimensioned at N. C C N4 C is the total number of observations or cells in the data C C NCLASS C is an integer expression which defines the number of classifications C for the analysis of variance where 1 <= NCLASS <= 7 C C NCELL C is an integer expression giving the number of observations C per cell in the analysis. C NCELL must be >= 1. C If NCLASS = 1, NCELL must be > 1 C C LEVEL C is an integer vector containing the number of levels C for each classification. C LEVEL must be dimensioned = NCLASS. C Each LEVEL(I) must be >= 2 C C VARNAM C is a vector of alphanumeric variables dimensioned >= NCLASS. C It should contain the names of the different factors C or classifications using up to 4 alphanumeric characters for each C C NPOOL SETS THE INTERaCTION LEVEL WHICH IS TO BE USED AS THE ERROR C TERM, ALL INTERACTIONS NPOOL AND ABOVE ARE POOLED AND C USED AS THE ERROR TERM TO GIVE A representative F STATISTIC C MEANS C is an integer vector which defines which sets of means C are to be printed. C The maximum number of main effects or interactions C for which means will be calculated is 10. C If fewer than 10 are required, the last number in the array should be 0. C The integer constants which define a main effect or interaction C are obtained as follows. C Suppose in a 5 way classification where the factors C are A,B,C,D, and E, the means for B, AC, and BDE are required. C Write the factor names in reverse order, place a 1 under each factor C required in the effect and a 0 under each preceeding factor C not required and convert the binary integer thus obtained C to a decimal integer. C C NOTE: C MEANS must be dimensioned at 10. C C E D C B A C 1 0 = 2 C 2 10 C C 1 0 1 = 5 C 2 10 C C 1 1 0 1 0 = 26 C 2 10 C C The vector MEANS could then be defined in the calling programme C as follows: C C DIMENSION MEANS(10) C DATA MEANS/2,5,26,7*0/ C C Interactions of higher order than 3 way are not permitted. C Improper values in the array MEANS are skipped. C C T and TTT C are floating point vectors required for work space. C They must be dimensioned to provide space greater than C or equal to the number of cells in the data. C C MDF C Stores the defrees of freedom and must be dimensioned C NCLASS C i.e. pi LEVEL(I) C I=1 C C If several analyses are being done, C T should be dimensioned in terms of the largest. C C REMARKS: C C 1. No output statements are required. C 2. Since the input data stored in the X vector is replaced C by residuals, the programme may not call ANOVA twice C unless the X vector is reinitialized. C 3. If MEANS(1) = 0, only the source of variation table is found. C C STORAGE OF DATA: C C The method of storing the data is in a singly-dimensioned array C where the subscripts corresponding to an observation must be C given an equivalent single subscript. C The subscripts of the first factor vary most rapidly C while those of the last factor vary least rapidly. C C EXAMPLE: C C Suppose you are investigating a design of three factors C (FA, FB, and FC) with factor FA at 3 levels C and factors FB and FC at 2 levels, C i.e. you have a 3 x 2 x 2 array which will be stored C linearly in a vector dimensioned at 12 with subscripts referring to C (a b c ,a b c ,a b c ,a b c ,a b c ,a b c ,a b c ,a b c ,a b c , C 1 1 1 2 1 1 3 1 1 1 2 1 2 2 1 3 2 1 1 1 2 2 1 2 3 1 2 C C a b c ,a b c ,a b c ) C 1 2 2 2 2 2 3 2 2 C C FA a a a C 1 2 3 C C FB b b b b b b C 1 2 1 2 1 2 C C FC c 4 2 5 6 2 4 C 1 C C c 6 2 3 1 4 2 C 2 C C Thus, this particular data would have to be input: C C 4 5 2 C 2 6 4 C 6 3 4 C 2 1 2 C C METHOD: C The means and anlysis of variance table are obtained with ANVAMI. C All calculations are done in single precision. C C TIMING: C Approximately N*(2**(NCLASS-2))/(NCELL*100) seconds where N is C the total number of observations. C C SOURCE: C Coded in FORTRAN IV by J. Ogilvie, University of Toronto, September, 1965. C Revised by Barbara Gibbins, Consulting Services, June, 1970. C C***************************************************************************** C C SAMPLE PROGRAM FOR ANOVA C C THIS PROGRAM USES THE TABLE OF DATA ILLUSTRATED UNDER THE C HEADING 'STORAGE OF DATA'. IF WE WISH TO CHECK ALL INTER- C ACTIONS, NOT THAT 'MEANS' WILL BE CALCULATED AS FOLLOWS: C C B A C 1 = 1 C 1 0 = 2 C 1 0 0 = 4 C 1 1 = 3 C 1 1 0 = 6 C 1 1 1 = 7 C SINCE THE NUMBER OF MEANS IS LESS THAN 10, THE LAST ONE C MUST EQUAL ZERO. C***************************************************************************** C 1 DIMENSION X(12),LEVEL(3),VNAME(3),MEANS(10),T(12) C 2 N=12; NCLASS=3; NCELL=1 C 5 PRINT 99 C 6 READ 1, (X(I),I=1,N) C 7 PRINT 11,X C 8 READ 2, (LEVEL(I),I=1,NCLASS) C 9 PRINT 22,LEVEL C 10 READ 3,(VNAME(I),I=1,NCLASS) C 11 PRINT 33,VNAME C 12 READ 2, (MEANS(I),I=1,7 C 13 PRINT 44, (MEANS(I),I=1,7 C 14 CALL ANOVA(X,N,NCLASS,NCELL,LEVEL,NNAME,MEANS,T) C 15 STOP C 16 99 FORMAT('0',T11,'INPUT DATA'/) C 17 1 FORMAT(3F10.5) C 18 11 FORMAT(' ',T11,'DATA MATRIX',(T21,3F9.1/) C 19 2 FORMAT(7I5) C 20 22 FORMAT('0',T11,'LEVELS',I13,2I9) C 21 3 FORMAT(3A4) C 22 33 FORMAT('0',T11,'FACTORS',A14,2A9) C 23 44 FORMAT('0',T11,'MEANS',7I6) C 24 END C C INPUT DATA C C DATA MATRIX 4.0 5.0 2.0 C 2.0 6.0 4.0 C 6.0 3.0 4.0 C 2.0 1.0 2.0 C C LEVELS 3 2 2 C C FACTORS FA FB FC C C MEANS 1 2 4 3 6 7 0 C C MEANS FOR FACTOR FA C 3.5000 3.7500 3.0000 C C MEANS FOR FACTOR FB C 4.0000 2.8333 C C MEANS FOR FACTOR FC C 3.8333 3.0000 C C C INTERACTION OF FA AND FB C C MEANS FOR FB AT LEVEL 1 OF FA C 5.0000 2.0000 C C MEANS FOR FB AT LEVEL 2 OF FA C 4.0000 3.5000 C C MEANS FOR FB AT LEVEL 3 OF FA C 3.0000 3.0000 C C C INTERACTION OF FB AND FC C C MEANS FOR FC AT LEVEL 1 OF FB C 3.6667 4.3333 C C MEANS FOR FC AT LEVEL 2 OF FB C 4.0000 1.6667 C C C TRIPLE INTERACTION OF FA, FB AND FC C C MEANS FOR FC AT LEVEL 1 OF FA AND AT LEVEL 1 OF FB C 4.0000 6.0000 C C MEANS FOR FC AT LEVEL 1 OF FA AND AT LEVEL 2 OF FB C 2.0000 2.0000 C C MEANS FOR FC AT LEVEL 2 OF FA AND AT LEVEL 1 OF FB C 5.0000 3.0000 C C MEANS FOR FC AT LEVEL 2 OF FA AND AT LEVEL 2 OF FB C 6.0000 1.0000 C C MEANS FOR FC AT LEVEL 3 OF FA AND AT LEVEL 1 OF FB C 2.0000 4.0000 C C MEANS FOR FC AT LEVEL 3 OF FA AND AT LEVEL 2 OF FB C 4.0000 2.0000 C C C SOURCE OF VARIATION DF SUM OF SQUARES MEAN SQUARE EFFECT NUMBER C C FA 2 1.166667 0.583333 1 C FB 1 4.083333 4.083333 2 C FA FB 2 5.166667 2.583333 3 C FC 1 2.083333 2.083333 4 C FA FC 2 11.166667 5.583333 5 C FB FC 1 6.750000 6.750000 6 C FA FB FC 2 0.500000 0.250000 7 C C THE RELATIVE ERROR OF THE TOTAL OF THE SUMS OF SQUARES TO THE TOTAL C SUM OF SQUARES IS 0.204E-05 C SUBROUTINEANOVAP(X,N4,NCLASS,NCELL,LEVEL,VARNAM,MEANS,T,TTT,MDF,NP *OOL) DIMENSIONX(N4),LEVEL(NCLASS),VARNAM(NCLASS),MEANS(10),E(7),V(3),NQ *(7),T(N4),MASK(7) REAL TTT(1) INTEGER MDF(1) INTEGER N(2) LOGICAL TEST EQUIVALENCE(E(1),V(1)),(N1,NQ(1)),(N2,NQ(2)),(N3,NQ(3)),(V1,V(1)), *(V2,V(2)),(V3,V(3)) DATA BLANK/' '/,MASK/"1,"2,"4,"10,"20,"40,"100/ DATA N/' ','*'/ NT=0 IF(NPOOL.GT.NCLASS.OR.NPOOL.LT.1)WRITE(6,1013)NCLASS IF(NPOOL.GT.NCLASS.OR.NPOOL.LT.1)NPOOL=NCLASS IF(NCLASS.LT.1)GOTO501 IF(NCLASS.GT.7)GOTO502 IF(NCELL.LT.1)NCELL=1 IF((NCLASS.EQ.1).AND.(NCELL.EQ.1))GOTO505 NOBS=1 DO5I=1,NCLASS LEVELI=LEVEL(I) IF(LEVELI.LT.2)GOTO503 5 NOBS=NOBS*LEVELI NN=NOBS*NCELL NUMBER=2**NCLASS-1 G=0.0 DO30I=1,NN 30 G=G+X(I) G=G/NN IF(MEANS(1).NE.0)WRITE(6,1000)G IF(NCELL.EQ.1)GOTO60 CSS=0.0 FNCELL=NCELL DO40I=1,NN 40 CSS=CSS+(X(I)-G)**2 DO55I=1,NOBS S=0.0 DO50L=I,NN,NOBS 50 S=S+X(L) N1=NOBS+I S1=S/FNCELL DO53L=N1,NN,NOBS 53 X(L)=X(L)-S1 55 X(I)=S G=G*FNCELL 60 DO190I=1,10 NEFF=MEANS(I) IF(NEFF.LE.0)GOTO200 IF(NEFF.GT.NUMBER)GOTO190 K=0 DO1992N1N=1,7 IF((MASK(N1N).AND.NEFF).NE.0.)K=K+1 IF(K.GT.3)GOTO190 1992 CONTINUE CALLANVAMI(0,X,N4,NCLASS,NCELL,LEVEL,VARNAM,NEFF,T,NT,E,NDF,S) J=0 DO90L=1,K 70 J=J+1 IF((MASK(J).AND.NEFF).EQ.0.)GOTO70 80 NQ(L)=LEVEL(J) 90 V(L)=VARNAM(J) GOTO(100,120,140),K 100 WRITE(6,1001)V1 WRITE(6,1002)(T(M),M=1,NT) GOTO190 120 WRITE(6,1003)V1,V2 DO130L=1,N1 WRITE(6,1004)V2,L,V1 130 WRITE(6,1002)(T(M),M=L,NT,N1) GOTO190 140 WRITE(6,1005)V1,V2,V3 L=N1*N2 DO150I1=1,N1 DO150I2=1,N2 WRITE(6,1006)V3,I1,V1,I2,V2 J=I1+N1*(I2-1) 150 WRITE(6,1002)(T(M),M=J,NT,L) 190 CONTINUE 200 TSS=0.0 DO210I=1,NOBS XII=X(I)-G TSS=TSS+XII*XII 210 X(I)=XII TSS=TSS/NCELL N2=NOBS-1 G=TSS N1=NUMBER-1 IF(N1.EQ.0)GOTO221 SSS=0. NNDF=0 DO220I=1,N1 CALLANVAMI(-1,X,N4,NCLASS,NCELL,LEVEL,VARNAM,I,T,NT,E,NDF,S) N2=N2-NDF S1=S/NDF TTT(I)=S MDF(I)=NDF IF(NPOOL.EQ.0)GO TO 219 NFIX=1 NPOOLE=0 ITEST=I 300 NPOOLE=NPOOLE+MOD(ITEST,2) ITEST=ITEST/2 IF(ITEST.NE.0)GO TO 300 IF(NPOOLE.LT.NPOOL)GO TO 219 NFIX=2 NNDF=NNDF+NDF SSS=SSS+S 219 CONTINUE 220 G=G-S 221 S=0.0 DO215I=1,NOBS 215 S=S+X(I)**2 S = S/NCELL S1=S/N2 G=G-S NFIX=1 NPOOLE=0 ITEST=NUMBER 603 NPOOLE=NPOOLE+MOD(ITEST,2) ITEST=ITEST/2 IF(ITEST.NE.0)GO TO 603 IF(NPOOLE.GE.NPOOL)NFIX=2 MDF(NUMBER)=N2 TTT(NUMBER)=S NNDF=NNDF+N2 SSS=SSS+S SES=SSS/NNDF WRITE(6,1015) 1015 FORMAT('1SOURSE OF VARIATION ',36X,'DF',4X,'SUM OF SQUARES ',8X,'M 1EAN SQUARE',11X,'MS/POOLED',3X,'EFFECT NUMBER',/) DO 604 I=1,NUMBER DO 605 J=1,7 605 E(J)=BLANK ITEST=I NFIX=1 NIT=0 NPOOLE=0 606 NMOD=0+MOD(ITEST,2) NIT=NIT+1 IF(NMOD.EQ.1)E(NIT)=VARNAM(NIT) NPOOLE=NPOOLE+NMOD ITEST=ITEST/2 IF(ITEST.NE.0)GO TO 606 IF(NPOOLE.GE.NPOOL)NFIX=2 SEX=TTT(I)/MDF(I) F=SEX/SES 604 WRITE(6,1014)N(NFIX),E,MDF(I),TTT(I),SEX,F,I WRITE(6,1012)NPOOL,NNDF,SSS,SES 1014 FORMAT(' ',A1,7A7,T52,I8,F18.6,2XF18.6,2XF18.6,6XI4) IF(NCELL.EQ.1)GOTO230 S=CSS-TSS NDF=NOBS*(NCELL-1) S1=S/NDF WRITE(6,1009)NDF,S,S1 DO225I=1,NOBS S=X(I) S1=0.0 N1=NOBS+I DO222L=N1,NN,NOBS X(L)=X(L)-S 222 S1=S1+X(L) 225 X(I)=-S1 230 G=ABS(G/TSS) IF(G.GT.1.0E-7)WRITE(6,1010)G RETURN 1000 FORMAT('1GRAND MEAN=',F11.4) 1001 FORMAT('0',/,'0MEANS FOR FACTOR',1X,A6) 1002 FORMAT(5X10F12.4) 1003 FORMAT('0'/'0INTERACTION OF ',1XA6,' AND ',A6) 0123000 C ADD COMMA'S NEXT TWO CARDS COLUMN 36 26 1004 FORMAT('0MEANS FOR',1XA6,' AT LEVEL',I4,' OF ',A6) 0124000 1005 FORMAT('0'/'0TRIPLE INTERACTION OF ',A6,', ',A6,' AND ',A6) 1006 FORMAT('0MEANS FOR ',A6,' AT LEVEL',I4,' OF ',A6,' AND AT LEVEL ' *,I4,' OF ',A6) 1009 FORMAT(30X'WITHIN CELLS',8XI8,F18.6,3XF18.6) 1010 FORMAT('0THE RELATIVE ERROR OF THE TOTAL OF THE SUMS OF SQUARES TO * THE TOTAL SUM OF SQUARES IS',E12.3,//) 1012 FORMAT(//,' ','* INTERACTIONS ',I1,' WAY AND ABOVE POOLED',T52,I 18,F18.6,3XF18.6,//) 1013 FORMAT(//,' ','ERROR]] NPOOL MUST BE GREATER THAN 1 BUT LESS THAN *NCLASS. NPOOL IS SET TO ',I1,' FOR THIS ANALYSIS',//) 500 WRITE(6,2500) STOP 501 WRITE(6,2501) GOTO500 502 WRITE(6,2502) GOTO500 503 WRITE(6,2503)VARNAM(I) GOTO500 505 WRITE(6,2505) GOTO500 2500 FORMAT(' ANOVA HAS NOT BEEN EXECUTED'/'1') 2501 FORMAT(' THE NUMBER OF CLASSIFICATIONS IS LESS THAN 1') 2502 FORMAT(' THE NUMBER OF CLASSIFICATIONS IS GREATER THAN 7') 2503 FORMAT(' FACTOR',A7,' HAS LESS THAN 2 LEVELS') 2505 FORMAT(' A ONE WAY CLASSIFICATION WITH ONE OBSN PER CELL IS A SING *LE VARIANCE') END SUBROUTINEANVAMI(INDEX,X,N4,NCLASS,NCELL,LEVEL,VARNAM,NEFF,TOTAL, *NTOTAL,ENAME,NDF,SSQTOT) DIMENSIONX(N4),LEVEL(NCLASS),VARNAM(NCLASS),ENAME(7),TOTAL(N4), *NL(8),MASK(7),KSUB(7),TEST(7) LOGICALTEST,FIRST DATA BLANK/' '/,MASK/"1,"2,"4,"10,"20,"40,"100/ C DATABLANK/' '/,MASK/Z1,Z2,Z4,Z8,Z10,Z20,Z40/ NL(1)=1 IF((NCLASS.LT.1).OR.(NCLASS.GT.7))GOTO300 IF((NEFF.GT.(2**NCLASS-1)).OR.(NEFF.LT.1))GOTO301 IF(NCELL.LT.1)NCELL=1 C OBTAIN TOTAL NUMBER OF OBSERVATIONS NN NN=1 DO100I=1,NCLASS IF(LEVEL(I).LT.2)GOTO302 NN=NN*LEVEL(I) 100 ENAME(I)=BLANK C COMPUTE NUMBER OF BITS REQUIRED IN NEFF IQ=1 DO101I=1,NCLASS IQ=IQ*2 II=I IF(IQ.GT.NEFF)GOTO102 101 CONTINUE 102 NT=II FIRST=.TRUE. C SET ARRAYS NL AND TEST FOR REQUIREDCLASSIFICATIONS NDF=1 NTOTAL=1 DO110I=1,NT IF((MASK(I).AND.NEFF).EQ.0.)GOTO105 TEST(I)=.TRUE. LEVELI=LEVEL(I) NL(I+1)=NL(I)*LEVELI NDF=NDF*(LEVELI-1) NTOTAL=NTOTAL*LEVELI ENAME(I)=VARNAM(I) GOTO110 105 TEST(I)=.FALSE. NL(I+1)=NL(I) 110 CONTINUE C OBTAIN DIVISOR FOR MEANS DIV=(NN*NCELL)/NTOTAL C INITIALIZE TOTALS AND SUBSCRIPTS FOR X DO120I=1,NTOTAL 120 TOTAL(I)=0.0 125 DO130I=1,NCLASS 130 KSUB(I)=0 C FROM THE SUBSCRIPTS KSUB(K) THE SUBSCRIPT FOR THE RELEVANT TOTAL C IS OBTAINED. X IS ADDED TO THE TOTAL AND KSUB(K) IS INCREASE BY 1 DO160I=1,NN KT=1 DO140J=1,NT 140 IF(TEST(J))KT=KT+NL(J)*KSUB(J) IF(.NOT.FIRST)GOTO141 TOTAL(KT)=TOTAL(KT)+X(I) GOTO145 141 X(I)=X(I)-TOTAL(KT)*DIV 145 DO150K=1,NCLASS KSUB(K)=KSUB(K)+1 IF(KSUB(K).LT.LEVEL(K))GOTO160 150 KSUB(K)=0 160 CONTINUE IF(.NOT.FIRST)RETURN SSQTOT=0.0 DO170I=1,NTOTAL 170 SSQTOT=SSQTOT+TOTAL(I)*TOTAL(I) SSQTOT=SSQTOT/DIV IF(INDEX.GE.1)RETURN DO180I=1,NTOTAL 180 TOTAL(I)=TOTAL(I)/DIV IF(INDEX.GE.0)RETURN C FOR INDEX=-1 THE MEANS ARE SUBTRACTED FROM THE OBSERVATIONS X FIRST=.FALSE. DIV=NCELL GOTO125 300 WRITE(6,400) 400 FORMAT(' THENUMBER OF CLASSIFICATIONS IS NOT BETWEEN 1 AND 7') GOTO500 301 WRITE(6,401)NEFF 401 FORMAT(' THE EFFECT NUMBER',I4,'IS NOT APPROPRIATE') GOTO500 302 WRITE(6,402)I 402 FORMAT(' THE NUMBER OF LEVELS OF FACTOR',I4,' IS LESS THAN 2') 500 WRITE(6,501) CALL CLOSE(6) STOP 501 FORMAT(' ANVAMI HAS NOT BEEN EXECUTED') END