C  CONTRAST
c
c  This program tests the hypothesis that survival (or recovery) rates
c  S(1),...S(n) are equal or specific contrasts of S.
c
c   Input: N = number of survival rates
c          S(i) = survival rates
c          VARCOV(i,j) = variance-covariance matrix
c          C(i) = contrast vector
c
c   Output: CHI = chi-square statistic
c           DF = degrees of freedom
c           PROB = probability level
c
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (IDIM=70)
      INTEGER*2 CH2FLG,CH3FLG,SEFLG,PARM,GRP(IDIM),NPG(IDIM)
      CHARACTER*40 FNAME,CCON,TITLE*80
      DIMENSION S(IDIM),VARCOV(IDIM,IDIM),C(IDIM,IDIM),
     ,   CT(IDIM,IDIM),CS(IDIM),CST(IDIM),CVC(IDIM,IDIM),
     ,   TMP1(IDIM,IDIM),TMP2(IDIM),CHI(1,1),
     ,   A(IDIM,IDIM),AT(IDIM,IDIM),SORIG(IDIM),VCORIG(IDIM,IDIM)
      LOGICAL IN1
      EQUIVALENCE (AT,CT),(A,CVC)
      DATA CCON/'CON COn Con cON cOn coN con'/,CH3FLG/0/
C
    5 PRINT *,'Program CONTRAST [8/1/90]'
      PRINT *,'   Tests hypotheses about survival or recovery rates.'
      PRINT *,'   Method: See Sauer & Williams (1989) Generalized'
      PRINT *,'           Procedures For Testing Hypotheses About'
      PRINT *,'           Survival Or Recovery Rates., J. Wildl.'
      PRINT *,'           Manage. 53:137-142.'
      PRINT *,' '
      PRINT *,'MODIFICATIONS:'
      PRINT *,' 8/1/90 : Fixed bug where user enters s(i),se(i) in'
      PRINT *,' 2nd analysis (covariances from 1st analysis were used.)'
      PRINT *,' '
      PRINT *,'   Results from this program will be saved in a file'
      PRINT *,'   called CONTRAST.OUT.  If the survival rates and '
      PRINT *,'   variance-covariance matrix are entered via the'
      PRINT *,'   keyboard, they will be saved in a file called'
      PRINT *,'   CONTRAST.TMP.  This file may be used as input to'
      PRINT *,'   CONTRAST in later runs.'
      PRINT *,' '
    6 IN1=.FALSE.
        DO 4 I=1,IDIM
          DO 4 J=1,IDIM
    4       VARCOV(I,J)=.0D0
      SEFLG=1
      PRINT *,'Enter input filename("CON" for keyboard input):'
      FNAME='                                        '
      READ (*,FMT='(A)')FNAME
      IF(FNAME(4:4).EQ.' '.AND.INDEX(CCON,FNAME(1:3)).GE.1)IN1=.TRUE.
      OPEN(1,FILE=FNAME,STATUS='OLD',ERR=500)
      IF(IN1)PRINT *,'Enter title to identify contrast:'
      IF(IN1)PRINT *,':'
      LINE=1
      READ(1,FMT='(A)',IOSTAT=IER,ERR=600)TITLE
    7 IF(IN1)PRINT *,'Enter number of survival rates(N):'
      LINE=2
      READ(1,*,IOSTAT=IER,ERR=600)N
      IF(N.GT.IDIM.OR.N.LE.1)THEN
        PRINT *,' '
        PRINT *,'Error: number of survival rates must be between 2 and '
     ,   ,IDIM
        GO TO 7
      ENDIF
      IF(IN1)PRINT 8
    8 FORMAT(/' Enter 1 if you will be entering standard errors,'/4X,
     ,'or 2 if you will be entering the variance-covariance matrix.'/)
      IF(IN1)PRINT *,'Choice (1 or 2):'
      LINE=3
      READ(1,*,IOSTAT=IER,ERR=600)SEFLG
    9 IF(SEFLG.LT.1.OR.SEFLG.GT.2)SEFLG=1
      IF(SEFLG.EQ.1)THEN
        DO 10 I=1,N
   11     IF(IN1)PRINT *,'Enter S,SE(',I,'):'
          LINE=LINE+1
          READ (1,*,IOSTAT=IER,ERR=600)S(I),SE
          IF(SE.LT..0D0)THEN
            PRINT *,'Error: Std. error less than zero, please retype-'
            GO TO 11
          ENDIF
   10     VARCOV(I,I)=SE*SE
      ELSE
        DO 15 I=1,N
          IF(IN1)PRINT *,'ENTER S(',I,'):'
          LINE=LINE+1
   15     READ (1,*,IOSTAT=IER,ERR=600)S(I)
        IF(IN1)PRINT 16
   16   FORMAT(/' Since the variance-covariance matrix is symmetric, onl
     1y the lower left half'/' needs to be entered.'/
     2'  ie.  var(S1)'/7x,'cov(S2,S1)  var(S2)'/
     3 7x,'cov(S3,S1)  cov(S3,S2)  var(S3) ...'/)
        DO 20 I=1,N
   19     IF(IN1)
     )    PRINT *,'Enter row ',I,' of the variance-covariance matrix (',
     ,     I,' values):'
          IF(IN1)PRINT *,'?'
          LINE=LINE+1
          READ (1,*,IOSTAT=IER,ERR=600)(VARCOV(J,I),J=1,I)
          IF(VARCOV(I,I).LT..0D0)THEN
            PRINT *,'Error: Variance less than zero, please retype-'
            GO TO 19
          ENDIF
   20     CONTINUE
      ENDIF
c                    ----------------------------------------------
c                      If keyboard input, save survival rates and
c                      variance-covariance matrix so they can be
c                      used for input later.
c                    ----------------------------------------------
      IF(IN1)THEN
        PRINT *,'Survival rates and variance-covariance matrix will be'
        PRINT *,'saved in a file called CONTRAST.TMP (which can be used'
        PRINT *,'as input to later runs of CONTRAST).'
        OPEN(3,FILE='CONTRAST.TMP')
        WRITE(3,'(A/I2/I2/(D16.8))')TITLE,N,2,(S(I),I=1,N)
        DO 24 I=1,N
   24     WRITE(3,'(5D16.8)')(VARCOV(J,I),J=1,I)
        CLOSE(3)
      ENDIF
c                        ------------------------------------------
c                         Write output to screen and output file
c                        ------------------------------------------
      OPEN(2,FILE='CONTRAST.OUT')
      WRITE(*,22)TITLE
      WRITE(2,22)TITLE
   22 FORMAT(/1X,A//' Null Hypothesis: Homogeneous survival rates'/)
      DO 21 I=1,N
        SORIG(I)=S(I)
        DO 21 J=1,I
          VARCOV(I,J)=VARCOV(J,I)
          VCORIG(I,J)=VARCOV(I,J)
   21     VCORIG(J,I)=VARCOV(I,J)
      NORIG=N
   40 WRITE(*,41)(S(I),I=1,N)
      WRITE(2,41)(S(I),I=1,N)
   41 FORMAT(/' Survival rates:'/(1X,8F9.4))
      IF(N.EQ.2)THEN
         WRITE(*,42)DABS(S(1)-S(2))
         WRITE(2,42)DABS(S(1)-S(2))
   42    FORMAT(/'  Difference between survival rates = ',F9.4)
      ENDIF
      WRITE(*,45)
      WRITE(2,45)
   45 FORMAT(/' ============== Variance-covariance matrix ===========')
      XN=N
      XN1=-.1D1/XN
      DO 55 I=1,N
        IF(ICHOIC.NE.2)THEN
          DO 50 J=1,N
   50       C(I,J)=XN1
          C(I,I)=.1D1+XN1
        ENDIF
        WRITE(*,56)(VARCOV(I,J),J=1,I)
   55   WRITE(2,56)(VARCOV(I,J),J=1,I)
   56 FORMAT(/(1X,8F9.6))
      NN1=N-1
c               Compute chi-square value:
c                                               -1
c                    = C * S * [C * VC * C'] * C * S
c
c               where C = contrast vector,
c                     S = vector of survival rates, and
c                    VC = variance-covariance matrix of survival rates.
c
   70 CALL MMULT(C,VARCOV,TMP1,NN1,N,N,IDIM,IDIM,IDIM)
      CALL TRANSP(C,CT,N,N)
      CALL MMULT(C,S,CS,NN1,N,1,IDIM,IDIM,1)
      CALL MMULT(TMP1,CT,CVC,NN1,N,NN1,IDIM,IDIM,IDIM)
      CALL MATINV(CVC,NN1)
      CALL MMULT(CS,CVC,TMP2,1,NN1,NN1,1,IDIM,IDIM)
      CALL MMULT(TMP2,CS,CHI,1,NN1,1,1,IDIM,1)
C
      WRITE(*,80)CHI,NN1,CHIPRB(NN1,CHI)
      WRITE(2,80)CHI,NN1,CHIPRB(NN1,CHI)
   80 FORMAT(/25X,'CHI-SQUARE VALUE   = ',F10.4/
     /        25X,'DEGREES OF FREEDOM = ',I5/
     /        25X,'PROBABILITY        = ',F10.4)
      DO 300 I=1,NORIG
        S(I)=SORIG(I)
        DO 300 J=1,NORIG
          VARCOV(J,I)=VCORIG(J,I)
  300     C(J,I)=.0D0
c
      N=NORIG
  304 IF(.NOT.IN1)THEN
        READ(1,*,END=399)ICHOIC
      ELSE
        WRITE(*,305)
  305   FORMAT(/' Would you like to '//
     1        '   1) compare survival rates with 2 or more groups'/
     2        '   2) enter contrast vectors to test a hypothesis'/
     3        '   3) start a new problem'/ '   4) quit'/)
        PRINT *,'Enter 1,2,3, or 4:'
        ICHOIC=99
        READ (*,*,END=400,ERR=306)ICHOIC
  306   IF(ICHOIC.LT.1.OR.ICHOIC.GT.4)THEN
          PRINT *,' -- Please enter a number between 1 and 4 --'
          GO TO 304
        ENDIF
      ENDIF
      IF(ICHOIC.EQ.4)GO TO 400
      IF(ICHOIC.EQ.3)THEN
        CLOSE(1)
        GO TO 5
      ENDIF
      IF(ICHOIC.EQ.2)THEN
  330   IF(IN1)THEN
          PRINT *,' '
          PRINT *,'Enter number of vectors:'
          READ(*,*,END=400,ERR=331,IOSTAT=IER)NN1
          DO 1330 I=1,NN1
            PRINT *,'Enter contrast vector # ',I,':'
 1330       READ(*,*,END=400,ERR=331,IOSTAT=IER)(C(I,J),J=1,N)
        ELSE
          READ(1,*)NN1
          DO 1331 I=1,NN1
 1331       READ(1,*,END=399)(C(1,J),J=1,N)
        ENDIF
        GO TO 332
  331   PRINT *,'Error reading contrast vector, please retype.(',IER,')'
        GO TO 330
  332   WRITE(*,*)' Input contrast vector:'
        WRITE(2,*)' Input contrast vector:'
        DO 1332 I=1,NN1
          WRITE(*,'(1X,10F7.3)')(C(I,J),J=1,N)
 1332     WRITE(2,'(1X,10F7.3)')(C(I,J),J=1,N)
      ELSE
        IF(CH3FLG.EQ.0)PRINT 365
  365 FORMAT(/' To compare 2 or more groups of survival rates, you must
     1specify which '/' group each survival rate belongs to. When prompt
     2ed for group number, enter'/' 0 for survival rates not in the cont
     3rast, 1 for survival rates in the first'/' group, 2 for survival r
     4ates in the 2nd group, ...etc.')
        CH3FLG=1
  370   IF(IN1)THEN
           PRINT *,' '
           PRINT *,'Enter a group number for each survival rate(',N,
     ,      ' values)'
           PRINT *,'?'
           READ(*,*,ERR=371,IOSTAT=IER)(GRP(I),I=1,N)
        ELSE
           READ(1,*,ERR=371,IOSTAT=IER)(GRP(I),I=1,N)
        ENDIF
        GO TO 372
  371   PRINT *,'Error reading group numbers, please retype. (',IER,')'
        GO TO 370
  372   NGRPS=0
        WRITE(*,373)(GRP(I),I=1,N)
        WRITE(2,373)(GRP(I),I=1,N)
  373   FORMAT(/' Input group numbers:'/(2X,20I3))
        DO 375 I=1,N
          IF(GRP(I).LT.0.OR.GRP(I).GT.N)THEN
            PRINT *,'Error: Group number must be between 0 and ',N
            GO TO 370
          ENDIF
  375     NPG(I)=0
        DO 380 I=1,N
          IF(GRP(I).GT.NGRPS)NGRPS=GRP(I)
          IF(GRP(I).GT.0)NPG(GRP(I))=NPG(GRP(I))+1
          DO 380 J=1,N
  380       A(J,I)=.0D0
        IF(NGRPS.LE.1)THEN
          PRINT *,'Error: There must be at least two groups specified'
          GO TO 370
        ENDIF
C
        DO 390 I=1,N
  390     IF(GRP(I).GT.0)A(GRP(I),I)=.1D1/FLOAT(NPG(GRP(I)))
C
        CALL MMULT(A,SORIG,S,NGRPS,N,1,IDIM,IDIM,1)
        CALL MMULT(A,VCORIG,TMP1,NGRPS,N,N,IDIM,IDIM,IDIM)
        CALL TRANSP(A,AT,NGRPS,N)
        CALL MMULT(TMP1,AT,VARCOV,NGRPS,N,NGRPS,IDIM,IDIM,IDIM)
        N=NGRPS
        GO TO 40
      ENDIF
      GO TO 70
  399 IN1=.TRUE.
      GO TO 304
  400 STOP
  500 PRINT *,'Error: file not found'
      GO TO 6
  600 PRINT *,'Error reading line ',LINE,' from input file (',IER,')'
      CLOSE(1)
      STOP
      END
      SUBROUTINE TRANSP(X,XT,N1,N2)
      PARAMETER (IDIM=70)
      DOUBLE PRECISION X(IDIM,IDIM),XT(IDIM,IDIM)
      DO 10 I=1,N1
        DO 10 J=1,N2
   10     XT(J,I)=X(I,J)
      RETURN
      END
      SUBROUTINE MMULT(A,B,C,N1,N2,N3,IA,IB,IC)
      DOUBLE PRECISION SUM,A(IA),B(IB),C(IA)
      DO 10 I=1,N1
        DO 10 J=1,N3
          SUM=.0D0
          DO 5 K=1,N2
    5       SUM=SUM+A(I+IA*(K-1))*B(K+IB*(J-1))
   10     C(I+IA*(J-1))=SUM
      RETURN
      END
      SUBROUTINE MATINV(A,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (IDIM=70)
      DIMENSION A(IDIM,IDIM),IPIVOT(IDIM),PIVOT(IDIM),INDEX(IDIM,2)
C
C     SUBROUTINE MATINV FROM MULTIVARIATE PROCEDURES FOR THE
C     BEHAVIOURAL SCIENCES BY W.W. COOLEY AND R.R. LOHNES, WILEY,
C     NEW YORK, 1962.
C
C     PROGRAMMED BY BURTON S.GARBOW,ARGONNE NATIONAL ALBORATORY
C     AND REPORTED IN IBM 704-709 SHARE LIBRARY AS AN F402.
C     THIS SUBROUTINE COMPUTES THE INVERSE AND DETERMINANT OF
C     MATRIX A, ORDER N, BY THE GAUSS-JORDAN METHOD. A-INVERSE
C     REPLACES A, AND THE DETERMINANT OF A IS PLACED IN DETERM.
C
      IF(N.EQ.1)THEN
        A(1,1)=.1D1/A(1,1)
        RETURN
      ENDIF
C     INITIALIZATION
      DETERM=.1D1
      DO 20 I=1,N
   20   IPIVOT(I)=0
C
C     SEARCH FOR PIVOT ELEMENT
C
      DO 550 I=1,N
         AMAX=.0D0
         DO 105 J=1,N
   50       IF (IPIVOT(J).EQ.1) GO TO 105
            DO 100 K=1,N
               IF (IPIVOT(K)-1) 80,100,740
   80          IF (DABS(AMAX).GE.DABS(A(J,K))) GO TO 100
               IROW=J
               ICOLUM=K
               AMAX=A(J,K)
  100          CONTINUE
  105       CONTINUE
         IPIVOT(ICOLUM)=IPIVOT(ICOLUM)+1
C
C     INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL
C
         IF (IROW.EQ.ICOLUM) GO TO 260
         DETERM=-DETERM
         DO 200 L=1,N
            SWAP=A(IROW,L)
            A(IROW,L)=A(ICOLUM,L)
            A(ICOLUM,L)=SWAP
  200       CONTINUE
  260    INDEX(I,1)=IROW
         INDEX(I,2)=ICOLUM
         PIVOT(I)=A(ICOLUM,ICOLUM)
C
C     DIVIDE PIVOT ROW BY PIVOT ELEMENT
C
         A(ICOLUM,ICOLUM)=.1D1
         DO 350 L=1,N
  350       A(ICOLUM,L)=A(ICOLUM,L)/PIVOT(I)
C
C     REDUCE NON-PIVOT ROWS
C
         DO 550 L1=1,N
            IF (L1.EQ.ICOLUM) GO TO 550
            T=A(L1,ICOLUM)
            A(L1,ICOLUM)=.0D0
            DO 450 L=1,N
               A(L1,L)=A(L1,L)-A(ICOLUM,L)*T
  450          CONTINUE
  550       CONTINUE
C
C     INTERCHANGE COLUMNS
C
      DO 710 I=1,N
         L=N+1-I
         IF(INDEX(L,1).EQ.INDEX(L,2)) GO TO 710
         JROW=INDEX(L,1)
         JCOLUM=INDEX(L,2)
         DO 705 K=1,N
            SWAP=A(K,JROW)
            A(K,JROW)=A(K,JCOLUM)
            A(K,JCOLUM)=SWAP
  705       CONTINUE
  710    CONTINUE
  740 RETURN
      END
      FUNCTION CHIPRB(IDF,CHSQ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C   COMPUTES THE PROB. OF THE CHI SQUARE VALUE 'CHSQ' WITH
C   'DF' DEGREES OF FREEDOM.
C
C   SUBROUTINE NEEDED: FUNCTION ALGAMA
C
      A=.5*FLOAT(IDF)
      X=.5*CHSQ
      CHIPRB=1.
      IF(X .LE. 0.001 .OR. A .LE. 0.49) GO TO 60
   10 TERM=1.
      SUM=0.
      COFN=A
      IF(X.GE.A.AND.X.GE.13.)GO TO 30
C     CONVERGENT SERIES FOR X .LT. A OR .LT. 13.
      FACT=-A
   20 TEMP=SUM
      SUM=SUM+TERM
      COFN=COFN+1.
      TERM=TERM*X/COFN
      IF(SUM-TEMP)50,50,20
C     ASYMPTOTIC SERIES FOR X .GTE. A AND X .GTE. 13.
   30 CHIPRB=0.
      FACT=X
   40 TEMP=SUM
      SUM=SUM+TERM
      COFN=COFN-1.
      RATIO=COFN/X
      TERM=TERM*RATIO
      IF(SUM-TEMP)50,50,40
   50 ARGMNT=DLOG(SUM)-X+A*DLOG(X)-ALGAMA(A)
      IF(ARGMNT .LT. -175.)GO TO 60
      CHIPRB=CHIPRB+EXP(ARGMNT)/FACT
   60 RETURN
      END
      FUNCTION ALGAMA(A)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C   FUNCTION COMPUTES THE LOG OF THE GAMMA FUNCTION GIVEN
C   THAT THE ARGUMENT IS AN INTEGER OR A MULTIPLE OF 1/2.
C
      ALGAMA=0.0
      I=INT(A)
      IF(INT(A+.6).NE.I)GO TO 20
      IF(A.LE. 2.0)GO TO 40
      DO 10 J=2,I-1
   10      ALGAMA=ALGAMA+ALOG(FLOAT(J))
      GO TO 40
C                  A IS A MULTIPLE OF 1/2
   20 ALGAMA=0.5723649429
      IF(A.LE. 0.5)GO TO 40
      DO 30 J=1,I*2-1,2
   30      ALGAMA=ALGAMA+ALOG(FLOAT(J)/2.)
   40 RETURN
      END