C PARTAREA.FOR C ESTIMATION OF THE AREA UNDER ONE ROC CURVE OVER AN A PRIORI C RANGE OF FALSE-POSITIVE RATES, AS WELL AS THE AREA UNDER TWO C CURVES OVER AN A PRIORI RANGE OF FALSE-POSITIVE RATES C C Reference: McClish, 1989. Analyzing a portion of the ROC curve, C Med. Decis. Making; 9: 190-195. C C This program uses an interactive input format. C IMPLICIT REAL(A-Z) DOUBLE PRECISION PARTA,CUM,PD,CCC1,CCC2,DERA,DERB + CON1,CON2,CON3,VARPA,CPA1PA2 LOGICAL FIRST CHARACTER*1 ONE,DEPEND,RERUN,ONLY PI=3.1415927 1 CONTINUE FIRST=.FALSE. C C INPUT THE DATA C WRITE(6,5) 5 FORMAT(' ONE ROC CURVE OR TWO? [1,2]') READ(5,6)ONE 6 FORMAT(A) WRITE(6,7) 7 FORMAT(' FROM THE MLE PROGRAMS, TYPE IN ESTIMATES OF') WRITE(6,8) 8 FORMAT(' A,VAR(A),B,VAR(B),COV(A,B) with a return betw each') READ *,A1,VARA1,B1,VARB1,COVA1B1 IF(ONE .EQ. '1')GOTO 20 WRITE(6,10) 10 FORMAT(' FOR THE SECOND ROC CURVE, TYPE IN ESTIMATES OF') WRITE(6,11) 11 FORMAT(' A2,VAR(A2),B2,VAR(B2),COV(A2,B2)') READ *,A2,VARA2,B2,VARB2,COVA2B2 WRITE(6,12) 12 FORMAT(' ARE THE TWO CURVES DEPENDENT? [Y,N]') READ(5,6)DEPEND IF(DEPEND .EQ. 'N' .OR. DEPEND .EQ. 'n')GOTO 20 WRITE(6,13) 13 FORMAT(' FROM THE MLE PROGRAMS, TYPE IN ESTIMATES OF') WRITE(6,14) 14 FORMAT( ' COV(A1,A2),COV(A1,B2),COV(B1,A2),COV(B1,B2)') 15 READ *,COVA1A2,COVA1B2,COVB1A2,COVB1B2 20 CONTINUE WRITE(6,21) 21 FORMAT(' TYPE THE RANGE OF FALSE-POSITIVE RATES') READ *,FP1,FP2 C C COMPUTE THE CUTOFF DECISION POINTS C C1=ZSCORE(FP1) C2=ZSCORE(FP2) IF(FP1 .LT. 0.5)C1=-C1 IF(FP1 .LE. 0.0)C1=-6.0 IF(FP2 .LT. 0.5)C2=-C2 IF(FP2 .GE. 1.0)C2=6.0 c print *,'c1,c2',c1,c2 C C COMPUTE THE PARTIAL AREA BY NUMERICAL INTEGRATION C AA=A1 BB=B1 VARA=VARA1 VARB=VARB1 COVAB=COVA1B1 25 CONTINUE PARTA=0.0 DO 50,N=C1,C2,0.0005 CUM=1.0-PNORM(AA+BB*N) PD=(1.0/SQRT(2.0*PI))*(EXP(-0.5*(N*N))) PARTA=PARTA+(CUM*PD*0.0005) 50 CONTINUE C C COMPUTE THE VARIANCE OF THE PARTIAL AREA C CC1=(C1+AA*BB/(1.0+BB*BB))*(SQRT(1.0+BB*BB)) CC2=(C2+AA*BB/(1.0+BB*BB))*(SQRT(1.0+BB*BB)) CCC1=(CC1*CC1)/2.0 CCC2=(CC2*CC2)/2.0 c print *,'c1p,c2p,c1pp,c2pp',cc1,cc2,ccc1,ccc2 CON1=EXP((-1.0*AA*AA)/(2.0*(1.0+BB*BB))) CON2=2.0*PI*(1.0+BB*BB) CON3=(1.0-PNORM(CC2))-(1.0-PNORM(CC1)) DERA=(CON1/SQRT(CON2))*CON3 DERB=((CON1/CON2)*(EXP(-1.0*CCC1)-EXP(-1.0*CCC2)))- + (((AA*BB*CON1)/(SQRT(2.0*PI)*SQRT((1.0+BB*BB)**3))) + *CON3) print *,DERA,DERB VARPA=DERA*DERA*VARA+DERB*DERB*VARB+2.0*DERA* + DERB*COVAB IF(ONE .EQ. '1')THEN C C OUTPUT FOR ONE ROC CURVE C WRITE(6,55)PARTA,VARPA 55 FORMAT(' THE PARTIAL AREA IS ',F6.4,' WITH VARIANCE ', + F10.8) GOTO 100 ENDIF C C IF TWO ROC CURVES, ESTIMATE AREA OF SECOND C IF(FIRST)GOTO 60 FIRST=.TRUE. PARTA1=PARTA VARPA1=VARPA DERA1=DERA DERB1=DERB AA=A2 BB=B2 VARA=VARA2 VARB=VARB2 COVAB=COVA2B2 GOTO 25 60 CONTINUE C C OUTPUT FOR TWO ROC CURVES C WRITE(6,65)PARTA1,VARPA1 65 FORMAT(' THE FIRST PARTIAL AREA IS ',F6.4,' WITH VARIANCE + ',F10.8) WRITE(6,66)PARTA,VARPA 66 FORMAT(' THE SECOND PARTIAL AREA IS ',F6.4,' WITH + VARIANCE ',F10.8) IF(DEPEND .EQ. 'N' .OR. DEPEND .EQ. 'n')GOTO 100 C C COMPUTE THE COVARIANCE BETW THE TWO DEPENDENT PARTIAL AREAS C CPA1PA2=DERA1*DERA*COVA1A2+DERA1*DERB*COVA1B2+DERB1* + DERA*COVB1A2+DERB1*DERB*COVB1B2 WRITE(6,70)CPA1PA2 70 FORMAT(' THE COVARIANCE BETWEEN THE TWO PARTIAL AREAS + IS ',F10.8) 100 CONTINUE WRITE(6,105) 105 FORMAT(' WOULD YOU LIKE TO RERUN THE PROGRAM? [Y,N]') READ(5,6)RERUN IF(RERUN .EQ. 'Y' .OR. RERUN .EQ. 'y')THEN WRITE(6,106) 106 FORMAT(' CHANGE THE FP RANGE ONLY? [Y,N]') READ(5,6)ONLY FIRST=.FALSE. IF(ONLY .EQ. 'Y' .OR. ONLY .EQ. 'y')GOTO 20 GOTO 1 ENDIF STOP END C*********************************************************************** C C FUNCTION ZSCORE (PVALUE) C ************************************************************ C C C WILL RETURN THE Z-SCORE OF A N(0,1) DISTRIBUTION C WHEN THE UPPER TAIL PROBABILITY IS GIVEN. C C C CHECKS THAT THE P-VALUE IS LESS THAN 1.0 AND GREATER THAN 0.0; C IF P-VALUE >= 1.0 IFAULT IS SET TO 1 (COMMON "BLKERR") C IF P-VALUE <= 0.0 IFAULT IS SET TO 2 (COMMON "BLKERR") C IF NO ERROR IFAULT IS ZERO C C THIS FUNCTION DETERMINES THE VALUE OF X SUCH THAT THE C C INTEGRAL OF EXP[ - (T**2)/2 ] C ----------------- C SQRT[ 2*PI ] C C EVALUATED FROM X TO INFINITY EQUALS Q. C WHERE Q IS THE CUMULATIVE PROBABILITY OF THE RIGHT C HAND TAIL. C C C THE FOLLOWING RATIONAL APPROXIMATION IS USED: C C C + C T + C T**2 C 0 1 2 C X = T - -------------------------- + ETA(Q) C 1 + D T + D T**2 + D T**3 C 1 2 3 C C WHERE THE ABSOLUTE VALUE OF ETAT(Q) < 4.5 X 10**-4 C AND T = SQRT [ LN (1/Q**2) ] C C0= 2.515517 ; D1= 1.432788 C C1= 0.802853 ; D2= 0.189269 C C2= 0.010328 ; D3= 0.001308 C C REF: ABRAMOWITZ AND STEGUN, HANDBOOK OF MATHEMATICAL C FUNCTIONS, NATIONAL BUREAU OF STANDARDS, 1968 C *********************************************************** C C IMPLICIT REAL (A-Z) INTEGER IFAULT COMMON /BLKERR/IFAULT C C0= 2.515517 D1= 1.432788 C1= 0.802853 D2= 0.189269 C2= 0.010328 D3= 0.001308 C SIGN= 1.0 ZSCORE= 0.0 IFAULT= 0 C... CHECKS THAT P-VALUE IS LESS THAN 1.0 AND GREATER THAN 0.0; C IF P-VALUE = 1.0 THEN IFAULT IS SET TO 1 IF(PVALUE.LT.1.0) GO TO 30 IFAULT= 1 RETURN C IF P-VALUE .LE. 0.0 THEN IFAULT IS SET TO 2 30 CONTINUE IF(PVALUE.GT.0) GO TO 40 IFAULT= 2 RETURN 40 CONTINUE IF(PVALUE.LE. 0.50) GO TO 50 SIGN= -1.0 PVALUE= 1.0 - PVALUE C 50 CONTINUE Q= PVALUE C T= SQRT( ALOG( 1/Q**2) ) C DEN= 1 + D1*T + D2*(T**2) + D3*(T**3) NUM= C0 + C1*T + C2*(T**2) C ZSCORE= (T - NUM/DEN)*SIGN C RETURN END C*********************************************************************** C C FUNCTION PNORM (ZSCORE) C ************************************************************************ C C C WILL RETURN THE UPPER TAIL PROBABILITY OF N(0,1) WHEN THE C CORRESPONDING ZSCORE IS GIVEN. C C C THIS FUNCTION DETERMINES THE VALUE OF Q SUCH THAT THE C C INTEGRAL OF EXP[ - (T**2)/2 ] C ----------------- C SQRT[ 2*PI ] C C EVALUATED FROM X TO INFINITY. C WHERE Q IS THE CUMULATIVE PROBABILITY OF THE RIGHT C HAND TAIL. C C C THE FOLLOWING RATIONAL APPROXIMATION IS USED: C C Q(X) = F(X) [ B T + B T**2 + B T**3 + B T**4 + B T**5 ] + ETA(X) C 1 2 3 4 5 C C WHERE THE ABSOLUTE VALUE OF ETAT(Q) < 7.5 X 10**-8 C AND T= 1/(1 + RX) ; R= 0.2316419 C C R= 0.2316419 C B1= 0.31938153 ; B2= -0.356563782 C B3= 1.781477937 ; B4= -1.821255978 C B5= 1.330274429 C F(X) IS THE INTEGRAND OF THE THE ABOVE INTEGRAL. C C REF: ABRAMOWITZ AND STEGUN, HANDBOOK OF MATHEMATICAL C FUNCTIONS, NATIONAL BUREAU OF STANDARDS, 1968 C ************************************************************************ C C IMPLICIT REAL (A-Z) C R= 0.2316419 PI= 3.14159268 B1= 0.31938153 B2= -0.356563782 B3= 1.781477937 B4= -1.821255978 B5= 1.330274429 C SIGN= 1.0 REF= 0.0 IF(ZSCORE.GE. 0.0) GO TO 50 SIGN= -1.0 REF= 1.0 ZSCORE= -1.0*ZSCORE C 50 CONTINUE X= ZSCORE C T= 1 / (1 + R*X) C POLY= B1*T + B2*T**2 + B3*T**3 + B4*T**4 + B5*T**5 FUNC= (1/SQRT(2*PI)) * EXP(-1.0*((X**2)/2)) C PNORM= SIGN*(FUNC * POLY) + REF C RETURN END