C GENERAL LINEAR HYPOTHESIS WITH CONTRASTS AUGUST 8, 1966 C THIS IS A SYSTEM/360 FORTRAN IV (LEVEL H) VERSION OF BMD06V C ORIGINALLY WRITTEN IN FORTRAN II. THE PROGRAM HAS BEEN SIFTED C AND SLIGHTLY MODIFIED TO MAKE IT OPERABLE ON THE 7094. C IT HAS BEEN FURTHER MODIFIED TO ALLOW ITS EXECUTION ON THE 360. C DIMENSION DATA(100),FMT(10),B(61),C(61,61),CI(61,61),NV(61), 1MV(61),NN(61),P(6),A(61,61),XM(61),GUNK(180),Q(6) DIMENSION DUMFMT(120) DOUBLE PRECISION P,RPR,GOBLE0,GOBLE1,GOBLE2,RRD,FMT, 1 RRE,URGLY,COEFFS COMMON DATA , GUNK , B , C , CI , NV COMMON MV , NN , Q , A , PR , RD COMMON NP , NO , NT , NEW , NH , NC COMMON ND , ON , RE , N1 , N2 , N3 COMMON N4 , NNT , NPP , MMT , MD , OON COMMON NTAPE , XM , NCOV , NT1 , NT2 C 902 FORMAT('1BMD06V - GENERAL LINEAR HYPOTHESIS WITH CONTRASTS - ', 1'REVISED SEPTEMBER 18, 1968'/ 241H HEALTH SCIENCES COMPUTING FACILITY, UCLA) C DATA P /6H66F1.0,6H33F2.0,6H22F3.0,6H16F4.0,6H13F5.0, 16H11F6.0/ NTAPE=5 CALL USAGEB('BMD06V') NT1=2 NT2=3 REWIND NT1 REWIND NT2 10 READ (5,900) RPR,RRD,NP,NO,NT,NEW,NH,NC,MD,MTAPE,NCOV,ND DATA GOBLE0/6HPROBLM/ IF(RPR.EQ. GOBLE0) GO TO 30 DATA GOBLE1/6HFINISH/ 15 IF(RPR.EQ. GOBLE1) GO TO 1000 WRITE(6,10000) 10000 FORMAT(87H0PROBLM OR FINISH IS PUNCHED WRONG ON THE PROPER CARD OR 1 EITHER CARD IS OUT OF SEQUENCE) STOP 6000 WRITE(6,7000) 7000 FORMAT(55H0COVARS IS PUNCHED WRONG OR COVARS CARD OUT OF SEQUENCE) STOP 8000 WRITE(6,9000) 9000 FORMAT(55H0COEFFS IS PUNCHED WRONG OR COEFFS CARD OUT OF SEQUEBCE) STOP 20 WRITE (6,901) STOP C C C 22 WRITE(6,1022) 1022 FORMAT(55H0VIOLATING THE LIMITATION FOR TOTAL NUMBER OF VARIABLES) STOP 66 WRITE(6,666) 666 FORMAT(41H0WRONG SET UP IN COL.13-72 ON CNTRST CARD) GO TO 170 2000 WRITE(6,3000) 3000 FORMAT(39H0ILLEGAL NUMBER OF CASES ON PRMBLM CARD) STOP 4000 WRITE(6,5000) 5000 FORMAT(45H0ILLEGALLY SPECIFIED NUMBER OF CNTRST CARD(S)) STOP 30 IF((MTAPE-NT1)*(MTAPE-NT2))8,20,8 8 CALL TPWD(MTAPE,NTAPE) 34 NNT=NP+NEW MMT=NNT-1 ON=NO IF(NNT.LE.0.OR.NNT.GT.60)GO TO 22 C WRITE(6,902) IF(NO.LE.0) GO TO 2000 IF(NH.LE.0) GO TO 4000 WRITE(6,903)RRD,NNT,NO IF(NC) 40, 40, 35 35 NPP=NP+1 GO TO 45 40 NPP=NP 45 DO 47 I=1,NNT 47 XM(I)=0.0 CALL READ IF(NEW) 50, 50, 1000 50 IF(NCOV) 53, 53, 52 52 READ(5,924)RRD,(NN(J), J=1,NCOV) DATA URGLY/6HCOVARS/ DATA COEFFS/6HCOEFFS/ IF(RRD.NE.URGLY) GO TO 6000 53 CALL ADJUST DO 170 IJ=1,NH READ(5,904)RRD,NT,ND,MD,(NV(I),I=1,MMT) DATA GOBLE2/6HCNTRST/ IF(RRD.EQ. GOBLE2) GO TO 60 55 WRITE (6,905) IJ GO TO 170 60 WRITE (6,906)IJ,(NV(I),I=1,MMT) NP=0 DO 70 I=1,MMT IF(NV(I)) 70, 70, 65 65 NP=NP+1 MV(NP)=I 70 CONTINUE IF(NP) 66, 66, 75 75 DO 80 J=1,NP N1=MV(J) DATA(J)=C(N1,NNT) DO 80 I=1,NP N2=MV(I) 80 CI(I,J)=C(N2,N1) WRITE (6,907) DO 85 I=1,NP WRITE (6,908)I 85 WRITE (6,909)(CI(I,J),J=1,NP),DATA(I) CALL INVERT (CI,NP,RD,NV,NN) WRITE (6,910) DO 90 I=1,NP WRITE (6,908)I 90 WRITE (6,909)(CI(I,J),J=1,NP) DO 95 J=1,NP B(J)=0.0 DO 95 I=1,NP 95 B(J)=B(J)+CI(I,J)*DATA(I) IF(MD) 98, 98, 96 96 WRITE (6,921) 98 RD=0.0 NDF=0 DO 108 I=1,NO READ(NTAPE)(DATA(J),J=1,NNT),PR RE=0.0 DO 100 J=1,NP N1=MV(J) 100 RE=RE+DATA(N1)*B(J) IF(RE) 101, 108, 101 101 OON=DATA(NNT)-RE NDF=NDF+1 IF(MD) 105, 105, 102 102 WRITE (6,922)I,DATA(NNT),RE,OON 105 RD=RD+(OON**2)*PR 108 CONTINUE REWIND NTAPE NDF=NDF-NP RE=NDF RD=RD/RE RE=SQRT(RD) DO 110 I=1,NP 110 DATA(I)=SQRT(RD*CI(I,I)) WRITE (6,911) WRITE (6,909)(B(I),I=1,NP) WRITE (6,912) WRITE (6,909)(DATA(I),I=1,NP) WRITE (6,913)RD WRITE (6,914)RE WRITE (6,923)NDF DATA FMT(1)/'( '/ DATA FMT(2)/'A6, '/ FMT(3)=P(ND) DATA FMT(4)/'/ '/ FMT(5) = FMT(1) DATA FMT(6)/'6X, '/ FMT(7)=P(ND) DATA FMT(8)/', '/ DATA FMT(9)/') '/ FMT(10) = FMT(9) DO 115 I=1,NT 115 READ(5,FMT)RRE,(A(I,J),J=1,NP) IF(RRE.NE.COEFFS) GO TO 8000 DO 120 I=1,NT DATA(I)=0.0 DO 120 J=1,NP 120 DATA(I)=DATA(I)+A(I,J)*B(J) DO 130 I=1,NT DUMFMT(I)=0.0 DO 125 K=1,NP DO 125 J=1,NP 125 DUMFMT(I)=DUMFMT(I)+A(I,K)*A(I,J)*CI(J,K) 130 DUMFMT(I)=SQRT(DUMFMT(I)*RD) DO 132 I=1,NT 132 NN(I)=I WRITE (6,915) N1=1 N2=0 135 IF(NT-7) 140, 140, 145 140 N2=N2+NT GO TO 150 145 N2=N2+7 150 WRITE (6,920)(NN(J),J=N1,N2) DO 160 I=1,NP 160 WRITE (6,916)MV(I),(A(J,I),J=N1,N2) WRITE (6,917)(DATA(J),J=N1,N2) WRITE (6,918) (DUMFMT(J),J=N1,N2) WRITE (6,919) NT=NT-7 IF(NT) 170, 170, 165 165 N1=N1+7 GO TO 135 170 CONTINUE GO TO 10 900 FORMAT(A6,A2,I2,I4,7I2,42X,I2) 901 FORMAT(33H1ERROR ON ASSIGNMENT OF TAPE UNIT) 903 FORMAT(17H0PROBLEM NUMBER A2//17H NO. OF VARIABLESI4/15H SAMPLE S 1IZE I6//) 904 FORMAT(A6,3I2,60I1) 905 FORMAT(43H0CNTRST IS PUNCHED WRONG ON CNTRST CARD NO.,I4,31H OR TH 1E CARD IS OUT OF SEQUENCE) 906 FORMAT(1H0//9H0CONTRAST/9H CARD NO.I3,5X,60I1) 907 FORMAT(1H0//23H SUMS OF CROSS PRODUCTS) 908 FORMAT(4H0ROWI3) 909 FORMAT(8F15.5) 910 FORMAT(1H0//31H INVERTED CROSS PRODUCTS MATRIX) 911 FORMAT(1H0//24H REGRESSION COEFFICIENTS) 912 FORMAT(47H0STANDARD DEVIATIONS OF REGRESSION COEFFICIENTS) 913 FORMAT(1H0//21H VARIANCE OF ESTIMATEF15.5) 914 FORMAT(23H STD. ERROR OF ESTIMATEF13.5//) 915 FORMAT(24H0CONTRASTS AND ESTIMATES) 916 FORMAT(1H I5,3X,7F15.4) 917 FORMAT(10H0ESTIMATE 7F15.5) 918 FORMAT(9H0STANDARD/10H ERROR OF 7F15.5) 919 FORMAT(9H ESTIMATE//) 920 FORMAT(10H0PARAMETER7X,I4,6(11X,I4)) 921 FORMAT(1H0/1H019X,18HTABLE OF RESIDUALS/12H OBSERVATION5X,7HY VALU 1E9X,10HY ESTIMATE8X,8HRESIDUAL) 922 FORMAT(I7,2F18.5,F16.5) 923 FORMAT(19H DEGREES OF FREEDOM5X,I6///) 924 FORMAT(A6,33I2,/(6X,33I2)) 9876 FORMAT(20A4) 1000 IF(NTAPE-5)1001,1002,1001 1001 REWIND NTAPE 1002 STOP END C SUBROUTINE ADJUST FOR BMD06V MAY 10, 1966 SUBROUTINE ADJUST DIMENSION DATA(100),FMT(180),B(61),C(61,61),CI(61,61),NV(61), 1MV(61),NN(61),P(6),A(61,61),XM(61) COMMON DATA , FMT , B , C , CI , NV COMMON MV , NN , P , A , PR , RD COMMON NP , NO , NT , NEW , NH , NC COMMON ND , ON , RE , N1 , N2 , N3 COMMON N4 , NNT , NPP , MMT , MD , OON COMMON NTAPE , XM , NCOV , NT1 , NT2 C MMMT=MMT-NCOV DO 3 I=1,MMT 3 MV(I)=I IF(NCOV) 40, 40, 10 10 WRITE (6,901) I=0 DO 5 J=1,MMT DO 1 K=1,NCOV IF(NN(K)-J)1,5,1 1 CONTINUE I=I+1 MV(I)=J 5 CONTINUE DO 20 I=1,NCOV K=NN(I) DO 15 L=1,MMMT J=MV(L) 15 DATA(L)=C(J,K)/C(J,J) WRITE (6,902)K 20 WRITE (6,900)(DATA(J), J=1,MMMT) 30 WRITE (6,903) GO TO 45 40 WRITE (6,904) 45 DO 50 L=1,MMMT J=MV(L) 50 DATA(L)=C(J,NNT)/C(J,J) WRITE (6,900)(DATA(J), J=1,MMMT) IF(NCOV) 1000, 1000, 60 60 DO 70 J=1,NNT DO 70 I=1,NNT 70 C(I,J)=0.0 DO 80 I=1,NCOV K=NN(I) 80 XM(K)=XM(K)/ON DO 100 I=1,NO READ(NT1)(DATA(J),J=1,NNT),PR DO 85 J=1,NCOV K=NN(J) 85 DATA(K)=DATA(K)-XM(K) DO 90 J=1,NNT DO 90 L=J,NNT 90 C(J,L)=C(J,L)+DATA(J)*DATA(L)*PR 100 WRITE(NT2)(DATA(J),J=1,NNT),PR REWIND NT1 END FILE NT2 REWIND NT2 NTAPE = NT2 N1=NNT-1 DO 110 I=1,N1 N2=I+1 DO 110 J=N2,NNT 110 C(J,I)=C(I,J) 900 FORMAT(8F15.5) 901 FORMAT(1H0//36H ORIGINAL MEANS FOR DESIGN VARIABLES) 902 FORMAT(13H0VARIABLE NO.I3,13H (COVARIATE)) 903 FORMAT(19H0DEPENDENT VARIABLE) 904 FORMAT(1H0//11H CELL MEANS) 9876 FORMAT(20A4) 1000 RETURN END C SUBROUTINE INVERT FOR BMD06V MAY 10, 1966 SUBROUTINE INVERT (A,N,D,L,M) DIMENSION A(61,61),L(61),M(61) C D=1.0 DO80 K=1,N L(K)=K M(K)=K BIGA=A(K,K) DO20 I=K,N DO20 J=K,N IF(ABS(BIGA)-ABS(A(I,J))) 10,20,20 10 BIGA=A(I,J) L(K)=I M(K)=J 20 CONTINUE J=L(K) IF(L(K)-K) 35,35,25 25 DO30 I=1,N HOLD=-A(K,I) A(K,I)=A(J,I) 30 A(J,I)=HOLD 35 I=M(K) IF(M(K)-K) 45,45,37 37 DO40 J=1,N HOLD=-A(J,K) A(J,K)=A(J,I) 40 A(J,I)=HOLD 45 DO55 I=1,N 46 IF(I-K)50,55,50 50 A(I,K)=A(I,K)/(-A(K,K)) 55 CONTINUE DO65 I=1,N DO65 J=1,N 56 IF(I-K) 57,65,57 57 IF(J-K) 60,65,60 60 A(I,J)=A(I,K)*A(K,J)+A(I,J) 65 CONTINUE DO75 J=1,N 68 IF(J-K)70,75,70 70 A(K,J)=A(K,J)/A(K,K) 75 CONTINUE D=D*A(K,K) A(K,K)=1.0/A(K,K) 80 CONTINUE K=N 100 K=(K-1) IF(K) 150,150,103 103 I=L(K) IF(I-K) 120,120,105 105 DO110 J=1,N HOLD=A(J,K) A(J,K)=-A(J,I) 110 A(J,I)=HOLD 120 J=M(K) IF(J-K) 100,100,125 125 DO130 I=1,N HOLD=A(K,I) A(K,I)=-A(J,I) 130 A(J,I)=HOLD GO TO 100 150 RETURN END C SUBROUTINE READ FOR BMD06V MAY 10, 1966 SUBROUTINE READ DIMENSION DATA(100),FMT(180),B(61),C(61,61),CI(61,61),NV(61), 1MV(61),NN(61),P(6),A(61,61),XM(61) DOUBLE PRECISION GOBLE3,RRD COMMON DATA , FMT , B , C , CI , NV COMMON MV , NN , P , A , PR , RD COMMON NP , NO , NT , NEW , NH , NC COMMON ND , ON , RE , N1 , N2 , N3 COMMON N4 , NNT , NPP , MMT , MD , OON COMMON NTAPE , XM , NCOV , NT1 , NT2 ASN(XX)=ATAN(XX/SQRT(1.0-XX**2)) C NEW=0 IF(NT) 20, 20, 8 8 OON=ON+1.0 WRITE (6,900) WRITE (6,901) DATA GOBLE3/6HTRNGEN/ DO 18 I=1,NT READ(5,905)RRD,(A(I,J), J=1,4) IF(RRD.EQ. GOBLE3) GO TO 12 10 WRITE (6,906) 11 STOP 12 N1=A(I,1) N2=A(I,2) N3=A(I,3) N4=A(I,4) IF(N2-8) 13, 14, 14 13 WRITE (6,902)I,N1,N2,N3 GO TO 18 14 IF(N2-10) 15, 16, 16 15 WRITE (6,903)I,N1,N2,N3,A(I,4) GO TO 18 16 WRITE (6,902)I,N1,N2,N3,N4 18 CONTINUE 20 IF(ND.GT.0.AND.ND.LE.10)GO TO 25 ND=1 WRITE(6,4000) 25 ND=ND*18 WRITE(6,5000) 5000 FORMAT(24H0VARIABLE FORMAT CARD(S)) READ (5,912)(FMT(I), I=1,ND) WRITE(6,6000)(FMT(I),I=1,ND) 6000 FORMAT(1X,18A4) DO 70 J=1,NNT DO 70 I=1,NNT 70 C(I,J)=0.0 ITAPE=5 IF(NTAPE.GT.0)ITAPE=NTAPE DO 280 I=1,NO 75 READ (ITAPE,FMT)(DATA(J), J=1,NPP) 77 IF(NC) 81, 81, 82 81 PR=1.0 GO TO 84 82 PR=DATA(NPP) IF(PR) 84, 81, 84 84 IF(NT) 260, 260, 90 90 DO 255 J=1,NT N1=A(J,1) N2=A(J,2) N3=A(J,3) D3=DATA(N3) N4=A(J,4) GO TO (110,110,110,140,150,160,170,180,190,200,210,220,230,240),N2 110 IF(D3)113,111,115 111 IF(N2-2) 114, 112, 114 112 D1=1.0 GO TO 250 113 WRITE (6,904)N2,N3,I,D3 NEW=NEW+1 114 D1=0.0 GO TO 250 115 GO TO (118,120,130), N2 118 D1=SQRT(D3) GO TO 250 120 D1=SQRT(D3)+SQRT(D3+1.0) GO TO 250 130 D1=ALOG10(D3) GO TO 250 140 D1=EXP(D3) GO TO 250 150 IF(-D3)154,154,113 154 RD=SQRT(D3) IF(RD-1.0) 155, 156, 113 155 D1=ASN(RD) GO TO 250 156 D1=1.5708 GO TO 250 160 IF(-D3)161,161,113 161 RD=SQRT(D3/OON) RE=SQRT((D3+1.0)/OON) IF(RD-1.0) 162, 163, 113 162 RD=ASN(RD) GO TO 164 163 RD=1.5708 164 IF(RE-1.0) 165, 166, 113 165 RE=ASN(RE) GO TO 167 166 RE=1.5708 167 D1=RD+RE GO TO 250 170 IF(D3)175,113,175 175 D1=1.0/D3 GO TO 250 180 D1=D3+A(J,4) GO TO 250 190 D1=D3*A(J,4) GO TO 250 200 D1=D3**N4 GO TO 250 210 D1=D3+DATA(N4) GO TO 250 220 D1=D3-DATA(N4) GO TO 250 230 D1=D3*DATA(N4) GO TO 250 240 IF(DATA(N4)) 242, 241, 242 241 WRITE (6,904)N2,N4,I,DATA(N4) NEW=NEW+1 GO TO 255 242 D1=D3/DATA(N4) 250 DATA(N1)=D1 255 CONTINUE 260 DO 270 J=1,NNT XM(J)=XM(J)+DATA(J) DO 270 L=J,NNT 270 C(J,L)=C(J,L)+DATA(J)*DATA(L)*PR 280 WRITE(NT1)(DATA(J),J=1,NNT),PR END FILE NT1 REWIND NT1 IF(NEW) 285, 285, 295 285 N1=NNT-1 DO 290 I=1,N1 N2=I+1 DO 290 J=N2,NNT 290 C(J,I)=C(I,J) 295 IF(MD) 1000, 1000, 296 296 WRITE (6,908) IF(NC) 297, 297, 298 297 N3=NNT GO TO 299 298 N3=NNT+1 299 DO 350 I=1,NO READ(NT1)(DATA(J),J=1,NNT),PR IF(NC) 305, 305, 300 300 DATA(N3)=PR 305 N1=1 N2=0 310 IF((N3-N2)-8) 315, 315, 320 315 N2=N3 GO TO 325 320 N2=N2+8 325 IF(N1-1) 330, 330, 335 330 WRITE (6,909)I,(DATA(J), J=N1,N2) GO TO 340 335 WRITE (6,910)(DATA(J), J=N1,N2) 340 IF(N3-N2) 350, 350, 345 345 N1=N1+8 GO TO 310 350 CONTINUE REWIND NT1 900 FORMAT(15H0TRANSFORMATION4X,12HVARIABLE NO.6X,7HTYPE OF9X,8HVARIAB 1LE6X,12H2ND VARIABLE) 901 FORMAT(12H CARD NO.9X,8HASSIGNED5X,14HTRANSFORMATION4X,11HTRANS 1FORMED4X,12HOR CONSTANT) 902 FORMAT(1H 4X,I4,13X,I4,13X,I3,13X,I3,10X,I6) 903 FORMAT(1H 4X,I4,13X,I4,13X,I3,13X,I3,10X,F12.5) 904 FORMAT(1H0//27H TRANSFORMATION OF THE TYPEI3,24H CANNOT BE PERFOR 1MED ON/9H VARIABLEI3,6H CASEI5,15H. THE VALUE ISF12.5) 905 FORMAT(A6,F3.0,F2.0,F3.0,F6.0) 906 FORMAT(29H0ERROR ON TRANSFORMATION CARD) 907 FORMAT(12I6) 908 FORMAT(1H0//1H 10X,34HINPUT DATA AS READ FROM CARDS/9H CASE N 1O.) 909 FORMAT(I6,8F13.5) 910 FORMAT(6X,8F13.5) 912 FORMAT (18A4) 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF 1IED, ASSUMED TO BE 1.) 9876 FORMAT(20A4) 1000 NTAPE=NT1 RETURN END C SUBROUTINE TPWD FOR BMD06V MAY 10, 1966 SUBROUTINE TPWD(NT1,NT2) IF(NT1)40,10,12 10 NT1=5 12 IF(NT1-NT2)14,19,14 14 IF(NT2.EQ.5)GO TO 18 17 REWIND NT2 19 IF(NT1-5)18,24,18 18 IF(NT1-6)22,40,22 22 REWIND NT1 24 NT2=NT1 28 RETURN 40 WRITE (6,49) STOP 49 FORMAT(25H ERROR ON TAPE ASSIGNMENT) END