C                     ROC_CONTGOLD.FOR
C
C	Last updated January 2004
C
C       THIS PROGRAM ESTIMATES THE ACCURACY OF ONE DIAGNOSTIC 
C	TEST OR COMPARES THE ACCURACY OF TWO DIAGNOSTIC
C       TESTS USING A MEASURE OF ACCURACY ANALOGOUS TO THE 
C       ROC AREA WHEN THE GOLD STANDARD (GS), OR TRUTH, IS
C       A CONTINUOUS-SCALE VARIABLE.  THE METHOD 
C       USES A NONPARAMETRIC APPROACH.  THE USER NEEDS TO 
C	MODIFY THIS CODE WHERE YOU SEE  "****" 
C
C	CURRENTLY, THE PROGRAM IS DIMENSIONED FOR TWO
C	DIAGNOSTIC TESTS AND UP TO 500 PATIENTS/SUBJECTS 
        implicit integer(a-z)
        real vcomp(2,500),area(2),truth(500),test(2,500),cov(2,2)
        real sum,zz,sum1,sum2,effectn,cutoff(15),binarea(2,15,15)
c
C	****ENTER THE INPUT FILE NAME HERE****
C	(Note that the input file should contain two fields, one for
C	the truth value and one for the diagnostic test result.  If you
C	are comparing the accuracies of two diagnostic tests, then
C	a third field is needed for the test result of the second
C	diagnostic.)
        open(unit=1,file='in.dat',status='old')
c
C	****ENTER THE OUTPUT FILE NAME HERE****
        open(unit=2,file='out.dat',status='new')
c
        print *,'Enter the total number of patients'        
        read *,total
c
	print *,'How many diagnostic tests? one(1) or two(2)'
	read *,modal
c
        do 30, k=1,total
	if(modal .eq. 1)read(1,55)truth(k),test(1,k)
        if(modal .eq. 2)read(1,55)truth(k),test(1,k),test(2,k)
 30     continue
c
C	****ENTER THE FORMAT OF THE INPUT FILE HERE****
 55     format(t55,f4.0,t66,f3.0,t73,f4.1)
c
C       COMPUTING THE STRUCTURAL COMPONENTS
c
        do 100, p=1,modal
c
        do 90,r=1,total
        sum=0.0
        do 85,q=1,total
        if(q .eq. r)goto 85
c
        if(truth(r) .eq. truth(q))sum=sum+0.5
        if((truth(r) .lt. truth(q)) .and. (test(p,r) .gt. test(p,q)))
     + sum=sum+1.0
        if((truth(r) .gt. truth(q)) .and. (test(p,r) .lt. test(p,q)))
     + sum=sum+1.0
        if((truth(r) .ne. truth(q)) .and. (test(p,r) .eq. test(p,q)))
     + sum=sum+0.5
c
 85     continue
        vcomp(p,r)=(sum/real(total-1))
 90     continue
 100    continue
c
C       COMPUTE THE ESTIMATES OF THE AREAS
        do 110,p=1,modal
        sum=0.0
        do 105,l=1,total
        sum=sum+vcomp(p,l)
 105    continue
        area(p)=sum/real(total)
 110    continue
c
C       COMPUTE THE ESTIMATED COVARIANCE MATRIX
        sum1=0.0
        sum2=0.0
        sum=0.0
c
        do 150, i=1,total
        sum1=sum1+(vcomp(1,i)-area(1))**2
	if(modal .eq. 1)goto 150
        sum2=sum2+(vcomp(2,i)-area(2))**2
        sum=sum+(vcomp(1,i)-area(1))*(vcomp(2,i)-area(2))
 150    continue
c
        effectn=real(total)/2.0
c
	cov(1,2)=0.0
        cov(1,1)=sum1/(effectn*(effectn-1.0))
	if(modal .eq. 1)goto 151
        cov(2,2)=sum2/(effectn*(effectn-1.0))
        cov(1,2)=sum/(effectn*(effectn-1.0))
        cov(2,1)=cov(1,2)
c
 151	continue
c
C       OUTPUT THE ESTIMATES
        do 300,p=1,modal
        write(2,225)p,area(p)
 225    format('TEST ',i1,' ESTIMATED ACCURACY: ',F6.4)
 300    continue
        write(2,301)
 301    format('ESTIMATED COVARIANCE MATRIX')
        do 310,p=1,modal
        write(2,303)(cov(p,a),a=1,2)
 303    format('    ',2(f9.6,2x))
 310    continue
	if(modal .eq. 1)goto 325
        zz=(area(1)-area(2))/sqrt(cov(1,1)+cov(2,2)-2.0*cov(1,2))
        write(2,320)zz
 320    format('Z-TEST STATISTIC IS ',f8.4)
 325	continue
c
c       
C	IF DESIRED, THE PROGRAM WILL ESTIMATE
C	PAIRWISE ROC AREAS FOR CATEGORIES OF THE
C	GOLD STANDARD WHICH YOU DEFINE
        print *,'Do you want a breakdown of accuracy' 
	print *,'across the gold standard range?'
        print *,'yes(1) or no(2)'
        read*,breakd
c
        if(breakd .eq. 2)goto 700
c
        print *,'How many categories do you want?'
        read *,category
        print *,'Enter the cutoff values with a return after each'
        do 405, k=2,category
        read *,cutoff(k)
 405    continue
c
        cutoff(1)=-99999999
        cutoff(category+1)=99999999
        do 500, k=1,category
        do 499, j=(k+1),category
c
        do 490, p=1,modal
        count=0
        sum=0.0
c
        do 480,r=1,total
        do 475,q=1,total
        if(q .eq. r)goto 475
        if((truth(r) .gt. cutoff(k+1)) .or. (truth(r) .lt. cutoff(k)))
     + goto 475
        if((truth(q) .gt. cutoff(j+1)) .or. (truth(q) .lt. cutoff(j)))
     + goto 475
c
        if(test(p,r) .lt. test(p,q))sum=sum+1.0
        if(test(p,r) .eq. test(p,q))sum=sum+0.5
        count=count+1
c
 475    continue
 480    continue
c
        binarea(p,k,j)=sum/real(count)
 490    continue
 499    continue
 500    continue
c
        write(2,510)
 510    format('ANALYSIS OF ACCURACY ALONG THE GOLD STANDARD RANGE')
        do 580, p=1,modal
        do 570, k=1,category
        do 560, j=(k+1),category
        write(2,600)p,k,j,binarea(p,k,j)
 560    continue
 570    continue
 580    continue
 600    format('test ',i1,' category ',i2,' vs category ',i2,': ',f6.4)
c
 700    continue
        stop
        end