PROGRAM ESTMAT C ESTIMATION OF SURVIVAL AND RECOVERY RATES UNDER FOUR DIFFERENT C STOCHASTIC MODELS. PARAMETER ESTIMATION IS VIA MAXIMUM LIKELIHOOD C SLIM = NUMBER YEARS OF BANDING C JLIM = NUMBER YEARS OF RECOVERY C RCVMAT = MATRIX OF RECOVERIES C AVEC = VECTOR OF BANDINGS C HEAD = ALPHA VECTOR OF TITLE C *** IPRINT CODE STRUCTURE: C *** 0 OR BLANK DEFAULT C *** 1 PAGE OF DEFINITIONS C *** OPTION CODE STRUCTURE: C *** 0 OR BLANK DEFAULT C *** 1 SPECIAL OUTPUT FOR MODEL 1 (SLOPES, ETC.) C *** 2 MODEL 3 (ABREVIATED) FOR VERY POOR DATA SETS C *** 3 SKIP DIRECTLY TO MODEL 1 C *** 4 SKIP TO MODEL 2 VIA MODEL 1 C *** 5 SKIP TO SAMPLE, FULL OUTPUT C *** 6 SKIP TO SAMPLE, OUTPUT SUPPRESSED C *** 7 SKIP TO SAMPLE, THEN TO MODEL 0 AND MODEL 1 INTEGER SLIM,YEAR1,OPTION DOUBLE PRECISION AVEC,RCVMAT COMMON /MAINCM/ RCVMAT(21,21),AVEC(20) CHARACTER*80 HEAD COMMON/CHERRY/ OPTION CHARACTER*12 OUTFIL WRITE(0,3) 3 FORMAT(' Enter the file name for the output'/) READ (5,'(A)') OUTFIL OPEN(6,FILE=OUTFIL,STATUS='UNKNOWN') C ENTER TITLE OF REPORT 7 WRITE(0,8) 8 FORMAT(/' Enter the TITLE between two apostrophes'/ 1' (example:''TITLE CARD'') (Use Ctrl-D to stop program)'/) 10 READ (5,*,END=45) HEAD C ENTER LIMITS ON YEARS WRITE(0,11) 11 FORMAT(/' Enter the following information (separated by ', 1'commas):', 2/' The number of years of banding, number of years of recovery,', 4/' and first year of banding study.', 5/' Then enter 1 if you wish a page of explanation, 0 if not.', 6/' Finally enter either 0 (for all models), ', 7 '2 (for very poor data),' 8/' 3 (for models 1, 2, and 3), or 4 (for models 2 and 3).'/) READ (5,*) SLIM,JLIM,YEAR1,IPRINT,OPTION IF (OPTION.EQ.5.OR.OPTION.EQ.6) CALL SAMPLE (OPTION,HEAD,IPRINT) IF (OPTION.EQ.5.OR.OPTION.EQ.6) GO TO 10 IF (OPTION.EQ.7) CALL SAMPLE(OPTION,HEAD,IPRINT) IF (OPTION.EQ.7) GO TO 10 IF (IPRINT.EQ.1.AND.OPTION.NE.2) CALL DEFINE DO 25 I=1,20 AVEC(I)=0.0 DO 25 JJ=1,20 RCVMAT(I,JJ)=0.0 25 CONTINUE C ENTER MATRIX OF RECOVERIES WRITE(0,12) SLIM 12 FORMAT(/' Enter each line of the matrix of recoveries,', 1' separating the numbers' 2/' on each line with commas. The matrix ', 2'is entered left justified'/' and should have',I3,' lines.'/) DO 35 II=1,SLIM WRITE(0,'(1X,A,I3,'':'')') ' Enter line',II READ (5,*) (RCVMAT(II,JJ),JJ=II,JLIM) 35 CONTINUE C ENTER VECTOR OF BANDINGS WRITE(0,13) 13 FORMAT(' Enter the vector of bandings, separating the values', 1' with commas.'/) READ (5,*) (AVEC(I),I=1,SLIM) WRITE(0,'(1x,A,A)') 'Program executing, output in file ',OUTFIL IF (OPTION.EQ.3) GO TO 43 IF (OPTION.EQ.4) GO TO 43 IF (OPTION.EQ.2) CALL SIMPLE(JLIM,SLIM,RCVMAT,AVEC,YEAR1,HEAD) IF (OPTION.EQ.2) GO TO 7 IF (JLIM.GT.2) THEN CALL MDL0(SLIM,JLIM,RCVMAT,AVEC,HEAD,YEAR1) ENDIF 43 CALL MDL1(SLIM,JLIM,RCVMAT,AVEC,HEAD,YEAR1,IPRINT,OPTION) 45 CLOSE(6) STOP 'ESTIMATE normal execution' END SUBROUTINE SAMPLE (OPTION,HEAD,IPRINT) C *** ESTIMATION OF SAMPLE SIZES FOR BANDING STUDIES ON ADULT BIRDS C *** TECHNIQUES DISCUSSED IN CHAPTER 9 OF BROWNIE ET AL C *** K = NUMBER OF YEARS OF BANDING C *** PHI = REASONABLE GUESS OF SURVIVAL C *** F = REASONABLE GUESS OF RECOVERY RATE C *** CVPHI = SPECIFIED PRECISION ON PHI C *** IFLAG = FLAG TO DESIGNATE LAST DATA CARD C *** OPTION=5; FULL OUTPUT C *** OPTION=6; SUPPRESSED OUTPUT C *** OPTION=7; SUPPRESSED OUTPUT, THEN CALL MODEL 0 AND MODEL 1 C IMPLICIT DOUBLE PRECISION (A-G,O-Z) INTEGER OPTION COMMON /SAMCOM/ PRMAT(21,21),R(21),C(21),T(21),TMINC(21),AVEC(20) CHARACTER*80 HEAD C C READ INPUT DATA C 5 READ (5,*) K,PHI,F,CVPHI,IFLAG WRITE (6,9) HEAD 9 FORMAT ('0',20X,A) WRITE (6,11) K 11 FORMAT ('0',' YEARS OF BANDING = ',I2) WRITE (6,12) PHI 12 FORMAT (' ',' SURVIVAL RATE = ',F5.4) WRITE (6,13) F 13 FORMAT (' ',' RECOVERY RATE = ',F5.4) WRITE (6,14) CVPHI 14 FORMAT (' ',' DESIRED COEFFICIENT OF VARIATION = ',F3.2) C C CALCULATE CONFIDENCE INTERVALS ABOUT MEAN AND WRITE THEM OUT. C SE=CVPHI*PHI SHI=PHI+1.96*SE SLO=PHI-1.96*SE WRITE (6,146) SLO,SHI 146 FORMAT (' ',' 95% CONFIDENCE INTERVALS ABOUT MEAN SURVIVAL = ',F5. '4,' - ',F5.4) IYEAR=0 C C INITIALIZE ARRAYS TO ZERO. C DO 16 I=1,K R(I)=0.0 C(I)=0.0 T(I)=0.0 TMINC(I)=0.0 DO 15 J=1,K PRMAT(I,J)=0.0 15 CONTINUE 16 CONTINUE C C GENERATE PROBABILITY MATRIX C PRMAT(1,1)=F DO 20 I=2,K PRMAT(1,I)=PHI*PRMAT(1,I-1) PRMAT(I,I)=PRMAT(1,1) 20 CONTINUE DO 30 I=1,K-1 DO 25 J=2,K-1 IF (I.GE.J) GO TO 25 PRMAT(I+1,J+1)=PRMAT(I,J) 25 CONTINUE 30 CONTINUE IF (OPTION.EQ.6) GO TO 74 IF (OPTION.EQ.7) GO TO 74 C C WRITE OUT PROBABILITY MATRIX C WRITE (6,40) 40 FORMAT ('0',' YEAR ',10X,'PROBABILITY OF BAND RECOVERY IN YEAR J') WRITE (6,50) 50 FORMAT (' ','BANDED') WRITE (6,60) 60 FORMAT ('+','______',10X,'______________________________________') IF (K.GT.10) GO TO 71 DO 70 I=1,K IYEAR =IYEAR+1 WRITE (6,65) IYEAR,(PRMAT(I,J),J=1,K) 65 FORMAT ('0',I4,5X,20(F6.4,5X)) 70 CONTINUE GO TO 74 71 IF (K.GT.15) WRITE (6,715) 715 FORMAT ('0','THE NUMBER OF YEARS OF BANDING IS TOO LARGE TO ', 1'PERMIT THE PRINTING OF THE PROBABILITY MATRIX') IF (K.GT.15) GO TO 74 DO 73 I=1,K IYEAR=IYEAR+1 WRITE (6,72) IYEAR,(PRMAT(I,J),J=1,K) 72 FORMAT ('0',I4,4X,20(F6.4,2X)) 73 CONTINUE 74 CONTINUE C C COMPUTE INTERMEDIATE STATISTICS. FIRST COMPUTE R(I), THE ROW TOTA C DO 80 I=1,K-1 DO 75 J=1,K IF (I.GT.J) GO TO 75 R(I)=PRMAT(I,J)+R(I) 75 CONTINUE 80 CONTINUE C C NOW COMPUTE C(I), THE COLUMN TOTALS. C C(1)=F DO 90 I=2,K-1 DO 85 J=1,K-1 C(I)=C(I)+PRMAT(J,I) 85 IF (I-J.EQ.1) GO TO 90 90 C(I)=C(I)+C(1) C C NOW COMPUTE T(I) AND T(I)-C(I). C DO 100 I=1,K-1 T(I)=(R(I)*C(I))/F TMINC(I)=T(I)-C(I) 100 CONTINUE IF (OPTION.EQ.6) GO TO 135 IF (OPTION.EQ.7) GO TO 135 C C WRITE R(I),C(I),T(I), AND T(I)-C(I) C WRITE (6,105) 105 FORMAT ('0','INTERMEDIATE STATISTICS') WRITE (6,110) 110 FORMAT ('0',' I R(I) C(I) T(I) T(I)-C(I)') WRITE (6,120) 120 FORMAT ('+','__________________________________________________') DO 130 I=1,K-1 WRITE (6,125) I,R(I),C(I),T(I),TMINC(I) 125 FORMAT ('0',I4,2(5X,F6.4),2(3X,F8.4)) 130 CONTINUE C C NOW FIND THE SUMS OF THE RECIPROCALS OF THE LAST TWO COLUMNS. C 135 A=0.0 B=0.0 DO 140 I=1,K-1 A=1.0/TMINC(I)+A B=1.0/T(I)+B 140 CONTINUE D=1.0/R(1)+(1.0/F)-2.0 C C COMPUTE THE FUNCTION H AND THE ESTIMATED SAMPLE SIZE N, C AND THE MATRIX OF EXPECTED RECOVERIES. C H=(A-B+D)/(K-1)**2 SAMSZE=H/(CVPHI**2) N=SAMSZE IF (OPTION.EQ.6) GO TO 143 DO 142 I=1,K DO 141 J=1,K PRMAT(I,J)=PRMAT(I,J)*N 141 CONTINUE 142 CONTINUE C C WRITE A, B, D, H, AND SAMSZE. C 143 WRITE (6,145) 145 FORMAT ('0','SUMMARY STATISTICS') WRITE (6,150) A 150 FORMAT ('0',' A = ',F9.4) WRITE (6,160) B 160 FORMAT (' ',' B = ', F9.4) WRITE (6,170) D 170 FORMAT (' ',' C = ',F9.4) WRITE (6,175) 175 FORMAT ('0','FUNCTION H AND ESTIMATED SAMPLE SIZE N') WRITE (6,180) H 180 FORMAT ('0',' H = ',F9.4) WRITE (6,190) N 190 FORMAT (' ',' N = ',I9) IF (OPTION.EQ.6) GO TO 280 IF (OPTION.NE.7) GO TO 199 IYEAR=1 DO 195 I=1,K 195 AVEC(I)=SAMSZE K1=K CALL MDL0(K,K1,PRMAT,AVEC,HEAD,IYEAR) CALL MDL1(K,K1,PRMAT,AVEC,HEAD,IYEAR,IPRINT,OPTION) GO TO 280 199 WRITE (6,200) 200 FORMAT ('0',' YEAR ',10X,'EXPECTED MATRIX OF RECOVERIES') WRITE (6,210) 210 FORMAT (' ','BANDED') WRITE (6,220) 220 FORMAT ('+','______',10X,'_____________________________') IYEAR=0 IF (K.GT.10) GO TO 240 DO 230 I=1,K IYEAR=IYEAR+1 WRITE (6,225) IYEAR,(PRMAT(I,J),J=1,K) 225 FORMAT ('0',I4,5X,20(F6.0,5X)) 230 CONTINUE GO TO 280 240 CONTINUE DO 260 I=1,K IYEAR=IYEAR+1 WRITE (6,265) IYEAR,(PRMAT(I,J),J=1,K) 265 FORMAT ('0',I2,2X,20(F6.0,1X)) 260 CONTINUE 280 CONTINUE IF (IFLAG.EQ.0) GO TO 5 RETURN END SUBROUTINE DECIDE(THCHI) C *** DECIDE OUTPUTS THE LIKELIHOOD RATIO TESTS AND ASSOCIATED C *** INFORMATION FOR MODELS 1, 2 AND 3. IMPLICIT DOUBLE PRECISION (A-G,O-Z) INTEGER OPTION DOUBLE PRECISION THCHI DIMENSION THCHI(100) COMMON/SNEAK/TEST,TEST2,TEST3,IDF,IDF2,IDF3 COMMON/CHERRY/ OPTION WRITE (6,5) 5 FORMAT ('0',30X,'TESTS OF VARIOUS MODELS AND ASSUMPTIONS'//) WRITE (6,10) 10 FORMAT (1X,'(IN EACH CASE THE NULL HYPOTHESIS BEING TESTED IS THA' 1,'T THE SIMPLEST MODEL, THE ONE WITH THE FEWEST PARAMETERS, ', 2 'FITS THE DATA)'//) WRITE (6,15) 15 FORMAT (1X,'TEST OF THE MODEL ASSUMING CONSTANT SURVIVAL AND RECO' 1,'VERY RATES AGAINST THE MODEL') IF (OPTION.EQ.4)GO TO 41 WRITE (6,20) 20 FORMAT (' ','ASSUMING TIME-SPECIFIC SURVIVAL AND RECOVERY RATES') WRITE (6,25) 25 FORMAT (' ','(F,S VS F(1),...,F(K) AND S(1),...,S(K-1) MODEL) 1 -- MODEL 3 VS. MODEL 1'/) WRITE (6,30) TEST 30 FORMAT (' ','CHI-SQUARED VALUE = ',F10.2) WRITE (6,35) THCHI(IDF) 35 FORMAT (' ','THEORETICAL CHI-SQUARE VALUE AT THE 5% LEVEL = ',F10. 12) WRITE (6,40) IDF 40 FORMAT (' ','DEGREES OF FREEDOM = ',I5) CALL CHI1(IDF,TEST) WRITE (6,45) 45 FORMAT (///' ',' ') WRITE (6,15) 41 WRITE (6,50) 50 FORMAT (1X,'ASSUMING CONSTANT SURVIVAL BUT TIME-SPECIFIC RECOVERY' 1,' RATES') WRITE (6,55) 55 FORMAT (1X,'(F,S VS F(1),...,F(K),S MODEL) -- MODEL 3 V' 1,'S. MODEL 2'/) WRITE (6,30) TEST2 WRITE (6,35) THCHI(IDF2) WRITE (6,40) IDF2 CALL CHI1(IDF2,TEST2) IF (OPTION.EQ.4) GO TO 71 WRITE (6,45) WRITE (6,60) 60 FORMAT (1X,'TEST OF THE MODEL ASSUMING CONSTANT SURVIVAL BUT TIME' 1,'-SPECIFIC RECOVERY RATES AGAINST THE MODEL') WRITE (6,65) 65 FORMAT (' ','ASSUMING TIME-SPECIFIC SURVIVAL AND RECOVERY RATES') WRITE (6,70) 70 FORMAT (' ','(F(1),...,F(K),S VS F(1),...,F(K) AND S(1),...,S(K- 11) MODEL) -- MODEL 2 VS. MODEL 1'/) WRITE (6,30) TEST3 WRITE (6,35) THCHI(IDF3) WRITE (6,40) IDF3 CALL CHI1(IDF3,TEST3) 71 CONTINUE RETURN END SUBROUTINE DEFINE C *** DEFINE PRINTS A TABLE OF DEFINITIONS AND INTRODUCTION (OPTIONAL) WRITE (6,5) 5 FORMAT ('0','THIS COMPUTER OUTPUT PRESENTS A DETAILED ANALYSIS OF' 1,' FOUR GENERAL STOCHASTIC MODELS FOR THE ANALYSIS OF BANDING') WRITE (6,10) 10 FORMAT (1X,'EXPERIMENTS. EACH MODEL IS BASED ON SEVERAL ASSUMPTI' 1,'ONS. IN EACH CASE PARAMETERS ARE ESTIMATED USING THE THEORY') WRITE (6,15) 15 FORMAT (1X,'OF MAXIMUM LIKELIHOOD AND THE ASSUMPTIONS ARE TESTED' 1,' STATISTICALLY. THE SEQUENCE PROGRESSES FROM A GENERAL MODEL,') WRITE (6,20) 20 FORMAT (1X,'BASED ON RELATIVELY FEW ASSUMPTIONS, TO A MODEL WHICH' 1,' MAKES VERY SIMPLE, BUT QUITE RESTRICTIVE ASSUMPTIONS.'//) WRITE (6,25) 25 FORMAT (' ','DEFINITIONS AND NOTATION') WRITE (6,30) 30 FORMAT ('+',25('_')/) WRITE (6,35) 35 FORMAT (' ',' K THE NUMBER OF YEARS OF BANDING'/) WRITE (6,40) 40 FORMAT (' ',' L THE NUMBER OF YEARS OF RECOVERY'/) WRITE (6,45) 45 FORMAT (1X,' F(I) THE RECOVERY RATE IN YEAR I. SPECIFICA' 1,'LLY, THE PROBABILITY THAT A BANDED BIRD IS RECOVERED') WRITE (6,50) 50 FORMAT (1X,' AND REPORTED TO THE BIRD BANDING LABORA' 1,'TORY IN YEAR I, GIVEN THAT IT WAS ALIVE AT THE BEGINNING') WRITE (6,55) 55 FORMAT (' ',' OF YEAR I.') WRITE (6,60) 60 FORMAT (' ',' I=1,2,...,K') WRITE (6,65) 65 FORMAT ('+',30X,'EXCEPT MODEL 0, WHERE I=2,3,...,K OR 2,3,...,K- 11'/) WRITE (6,70) 70 FORMAT (1X,' F THE CONSTANT RECOVERY RATE (NOT AN AVER' 1,'AGE VALUE)'/) WRITE (6,75) 75 FORMAT (1X,' S(I) THE SURVIVAL RATE IN YEAR I. SPECIFICA' 1,'LLY, THE PROBABILITY OF SURVIVAL OF A BIRD IN YEAR I,') WRITE (6,80) 80 FORMAT (1X,' GIVEN THAT IT WAS ALIVE AT THE BEGINNIN' 1,'G OF YEAR I.') WRITE (6,85) 85 FORMAT (' ',' I=1,2,...,K-1'/) WRITE (6,90) 90 FORMAT (1X,' ESTIMATES OF SURVIVAL PERTAIN TO THE PE' 1,'RIOD FROM THE TIME OF BANDING IN YEAR I TO THE TIME OF ') WRITE (6,95) 95 FORMAT (1X,' BANDING IN YEAR I+1'/) WRITE (6,100) 100 FORMAT (1X,' S THE CONSTANT SURVIVAL RATE (NOT AN AVER' 1,'AGE VALUE)'/) C WRITE (6,105) C 105 FORMAT (1X,' Y(I) AN ESTIMATE OF NONHUNTING MORTALITY RAT' C 1,'E IN YEAR I.') C WRITE (6,110) C 110 FORMAT (1X,' THIS IS USUALLY NOT ESTIMABLE AND IS EM' C 1,'PLOYED ONLY IN A SPECIAL SECTION OF MODEL 1') C WRITE (6,75) WRITE (6,105) 105 FORMAT (1X,' C(I) COLUMN TOTAL OF THE RECOVERY MATRIX IN' 1,' YEAR I -- THE TOTAL RECOVERIES IN CALENDAR YEAR I') WRITE (6,110) 110 FORMAT (1X,' I=1,2,...,L'/) WRITE (6,115) 115 FORMAT (1X,' R(I) ROW TOTAL OF THE RECOVERY MATRIX IN YEA' 1,'R I -- THE TOTAL RECOVERIES FROM BANDING IN YEAR I') WRITE (6,120) 120 FORMAT (' ',' I=1,2,...,K'/) WRITE (6,125) 125 FORMAT (1X,' T(I) A BLOCK TOTAL: THE TOTAL RECOVERIES ', 1'IN AND AFTER YEAR I FROM ALL BIRDS BANDED PRIOR TO, AND IN, ', 2'YEAR I') WRITE (6,120) WRITE (6,130) 130 FORMAT (' ',' N(I) NUMBER BANDED IN YEAR I') WRITE (6,120) WRITE (6,135) 135 FORMAT (1X,' RHO(I) TOTAL RECOVERY RATE FROM BIRDS BANDED I' 1,'N YEAR I') WRITE (6,120) RETURN END SUBROUTINE INVERT(A,N) C MATRIX INVERSION BY GAUSS-JORDAN ELIMINATION IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(21,21), B(21), C(21), LZ(21) DO 5 J=1,N 5 LZ(J)=J DO 55 I=1,N K=I Y=A(I,I) L=I-1 LP=I+1 IF (N-LP) 25,10,10 10 DO 20 J=LP,N W=A(I,J) IF (DABS(W)-DABS(Y)) 20,20,15 15 K=J Y=W 20 CONTINUE 25 DO 30 J=1,N C(J)=A(J,K) A(J,K)=A(J,I) A(J,I)=-C(J)/Y A(I,J)=A(I,J)/Y 30 B(J)=A(I,J) A(I,I)=1.0/Y J=LZ(I) LZ(I)=LZ(K) LZ(K)=J DO 50 K=1,N IF (I-K) 35,50,35 35 DO 45 J=1,N IF (I-J) 40,45,40 40 A(K,J)=A(K,J)-B(J)*C(K) 45 CONTINUE 50 CONTINUE 55 CONTINUE DO 80 I=1,N IF (I-LZ(I)) 60,80,60 60 K=I+1 DO 75 J=K,N IF (I-LZ(J)) 75,65,75 65 M=LZ(I) LZ(I)=LZ(J) LZ(J)=M DO 70 L=1,N C(L)=A(I,L) A(I,L)=A(J,L) 70 A(J,L)=C(L) 75 CONTINUE 80 CONTINUE RETURN END