C ROC_ORDGOLD.FOR C C FORTRAN PROGRAM TO ESTIMATE THE DIAGNOSTIC C ACCURACY OF A TEST(S) WHEN THE GOLD STANDARD C IS ORDINAL-SCALE. THE METHOD USES A C NONPARAMETRIC APPROACH. THE USER NEEDS TO C MODIFY THIS CODE WHERE YOU SEE "****" c c This version of the program can either estimate C the accuracy of a single diagnostic test or C compare 2 tests' accuracies. c c Written and documented by Nancy Obuchowski. Last revised January 2004. c c The program is currently dimensioned for up to 15 truth categories c and 500 patients per truth category. c implicit integer(a-z) dimension nn(15),data(2,15,500) real sum,area(2,15,15),prev(15),loss(15,15),cov(2,15,15,15,15) real summary1(2),summary2(2),summary3(2),denom1,denom2 real sqsum,vv(2,15,15,500),first,second, + cov12(15,15,15,15),ses1(2),ses2(2),ses3(2) real var(2,15,15),vars1(2),vars2(2),vars3(2), + covars1(2),covars2(2),covars3(2),cov12s1,cov12s2,cov12s3 c c ****ENTER THE INPUT AND OUTPUT FILE NAMES HERE **** c The input file is in unit 1. It must be sorted by truth state. open(unit=1,file='in.dat',status='old') open(unit=2,file='out.dat',status='new') c c ****Enter the # of observations per truth category **** c Note that here we have set up for 4 truth categories. nn(1)=114 nn(2)=21 nn(3)=19 nn(4)=87 c c ****Enter the population frequencies here**** c they must sum to 1.0. prev(1)=0.473 prev(2)=0.087 prev(3)=0.079 prev(4)=0.361 c c ****Enter values for penalty function here**** c if you don't want to assign different penalties c for different errors in misclassification, then c then give penalties all values of 1.0. loss(1,2)=1.0 loss(1,3)=1.0 loss(2,3)=1.0 loss(1,4)=1.0 loss(2,4)=1.0 loss(3,4)=1.0 c c c ****Enter the # of truth CATEGORIES here**** ii=4 c do 5, k=1,ii do 4, l=1,ii do 3, t=1,ii do 2, u=1,ii cov(1,k,l,t,u)=0.0 cov(2,k,l,t,u)=0.0 2 continue 3 continue 4 continue 5 continue c C ****ENTER THE NUMBER OF DIAGNOSTIC TESTS**** C ONE (1) OR TWO (2) tesst=1 c c c Read in the data do 35, t=1,ii do 10, k=1,nn(t) if(tesst .eq. 1)read(1,36)truth,data(1,t,k) if(tesst .eq. 2)then read(1,36)truth,data(1,t,k),data(2,t,k) endif 10 continue 35 continue c c ****ENTER FORMAT STATEMENT FOR INPUT FILE**** 36 format(t60,i1,t69,i1,2x,i2) c c do 300, modal=1,tesst c c Compute binary estimates between state k and state l c do 100, k=1,ii do 90, l=(k+1),ii c c Make comparisons based on kernel c do 60, m=1,nn(k) sum=0.0 do 59, n=1,nn(l) if(data(modal,k,m) .lt. data(modal,l,n))sum=sum+1.0 if(data(modal,k,m) .eq. data(modal,l,n))sum=sum+0.5 59 continue c c compute the structural component for m-th patient in k vs. state l c vv(modal,k,l,m)=sum/real(nn(l)) c c 60 continue c c do 70, n=1,nn(l) sum=0.0 do 69, m=1,nn(k) if(data(modal,k,m) .lt. data(modal,l,n))sum=sum+1.0 if(data(modal,k,m) .eq. data(modal,l,n))sum=sum+0.5 69 continue c c compute the structural component for n-th patient in l vs. state k c vv(modal,l,k,n)=sum/real(nn(k)) c c 70 continue c c Compute the binary area estimate c sum=0.0 do 75, m=1,nn(k) sum=sum+vv(modal,k,l,m) 75 continue area(modal,k,l)=sum/real(nn(k)) area(modal,l,k)=area(modal,k,l) c c c compute the variance for this binary estimate c sqsum=0.0 do 78, m=1,nn(k) sqsum=sqsum+(vv(modal,k,l,m)-area(modal,k,l))**2 78 continue first=sqsum/real(nn(k)-1) c sqsum=0.0 do 79, n=1,nn(l) sqsum=sqsum+(vv(modal,l,k,n)-area(modal,k,l))**2 79 continue second=sqsum/real(nn(l)-1) c c var(modal,k,l)=first/real(nn(k))+second/real(nn(l)) c c c 90 continue 100 continue c c c Compute the covariance between estimated areas in same modality c c *****Same firsts do 140, k=1,(ii-2) do 130, l=(k+1),(ii-1) c do 125, t=(l+1),ii do 122, m=1,nn(k) cov(modal,k,l,k,t)=cov(modal,k,l,k,t)+(vv(modal,k,l,m)- + area(modal,k,l))*(vv(modal,k,t,m)-area(modal,k,t)) 122 continue c cov(modal,k,l,k,t)=(cov(modal,k,l,k,t)/real(nn(k)-1))/ + real(nn(k)) c 125 continue c 130 continue 140 continue c c *****Same lasts if(ii .lt. 3)goto 171 do 170, k=3,ii do 160, l=1,(k-2) c do 155, t=(l+1),(k-1) do 150, m=1,nn(k) cov(modal,l,k,t,k)=cov(modal,l,k,t,k)+(vv(modal,k,l,m)- + area(modal,l,k))*(vv(modal,k,t,m)-area(modal,t,k)) 150 continue c cov(modal,l,k,t,k)=(cov(modal,l,k,t,k)/real(nn(k)-1))/ + real(nn(k)) c 155 continue c 160 continue 170 continue c 171 continue c c *****Same middles do 190, k=2,(ii-1) do 185, l=1,(ii-2) if(l .ge. k)goto 185 c do 180, t=(k+1),ii do 175, m=1,nn(k) cov(modal,l,k,k,t)=cov(modal,l,k,k,t)+(vv(modal,k,l,m)- + area(modal,l,k))*(vv(modal,k,t,m)-area(modal,k,t)) 175 continue c cov(modal,l,k,k,t)=(cov(modal,l,k,k,t)/real(nn(k)-1))/ + real(nn(k)) c 180 continue c 185 continue 190 continue c c c c Compute the summary indices of accuracy and their variances c summary1(modal)=0.0 denom1=0.0 summary2(modal)=0.0 denom2=0.0 summary3(modal)=0.0 vars1(modal)=0.0 vars2(modal)=0.0 vars3(modal)=0.0 c do 210, o=1,ii do 205, s=(o+1),ii denom1=denom1+real(nn(o)*nn(s)) denom2=denom2+prev(s)*prev(o) 205 continue 210 continue c print *,denom1 c do 220, o=1,ii do 215, s=(o+1),ii summary1(modal)=summary1(modal)+(real(nn(o)*nn(s))/ + denom1)*area(modal,o,s) vars1(modal)=vars1(modal)+((real(nn(o)*nn(s))/ + denom1)**2)*var(modal,o,s) summary2(modal)=summary2(modal)+(prev(s)*prev(o)/ + denom2)*area(modal,o,s) vars2(modal)=vars2(modal)+((prev(s)*prev(o)/ + denom2)**2)*var(modal,o,s) summary3(modal)=summary3(modal)+(prev(o)*prev(s)/ + denom2)*loss(o,s)*(1.0-area(modal,o,s)) vars3(modal)=vars3(modal)+((prev(o)*prev(s)/ + denom2)**2)*(loss(o,s)**2)*var(modal,o,s) 215 continue 220 continue c summary3(modal)=1.0-summary3(modal) c c c Sum up the covariances and multiply by 2 c covars1(modal)=0.0 covars2(modal)=0.0 covars3(modal)=0.0 c do 240, o=1,ii do 239, s=1,ii do 238, t=1,ii do 237, u=1,ii covars1(modal)=covars1(modal)+(real(nn(o)*nn(s))/ + denom1)*(real(nn(t)*nn(u))/denom1)*cov(modal,o,s,t,u)*2 c covars2(modal)=covars2(modal)+(prev(s)*prev(o)/denom2) + *(prev(t)*prev(u)/denom2)*cov(modal,o,s,t,u)*2 c covars3(modal)=covars3(modal)+(prev(o)*prev(s)/denom2) + *(prev(t)*prev(u)/denom2)*loss(o,s)*loss(t,u)* + cov(modal,o,s,t,u)*2 c 237 continue 238 continue 239 continue 240 continue c c 300 continue c c if(tesst .eq. 1)then write(2,301) 301 format(8x,'Area', 3x,'Variance') do 325, o=1,(ii-1) do 324, s=(o+1),ii write(2,320)o,s,area(1,o,s),var(1,o,s) 320 format(i2,' v',i2,2x,f6.4,1x,f9.7) 324 continue 325 continue c ses1(1)=sqrt(vars1(1)+covars1(1)) ses2(1)=sqrt(vars2(1)+covars2(1)) ses3(1)=sqrt(vars3(1)+covars3(1)) c write(2,352) 352 format(16x,'SummaryArea', 2x,'SE') c write(2,655)summary1(1),ses1(1) write(2,656)summary2(1),ses2(1) write(2,657)summary3(1),ses3(1) c goto 1000 endif c c Here, we compute the covariance between the two modalities c c **first, the same two disease groups in both modalities c do 400, k=1,ii do 390, l=(k+1),ii c sqsum=0.0 do 378, m=1,nn(k) sqsum=sqsum+(vv(1,k,l,m)-area(1,k,l))*(vv(2,k,l,m)- + area(2,k,l)) 378 continue first=sqsum/real(nn(k)-1) c sqsum=0.0 do 379, n=1,nn(l) sqsum=sqsum+(vv(1,l,k,n)-area(1,k,l))*(vv(2,l,k,n)- + area(2,k,l)) 379 continue second=sqsum/real(nn(l)-1) c c cov12(k,l,k,l)=first/real(nn(k))+second/real(nn(l)) c 390 continue 400 continue c c **second, only one common disease group in the two modalities c c *****Same firsts do 440, k=1,(ii-2) do 430, l=(k+1),ii c do 425, t=(k+1),ii if(t .eq. l)goto 425 do 422, m=1,nn(k) cov12(k,l,k,t)=cov12(k,l,k,t)+(vv(1,k,l,m)- + area(1,k,l))*(vv(2,k,t,m)-area(2,k,t)) 422 continue c cov12(k,l,k,t)=(cov12(k,l,k,t)/real(nn(k)-1))/ + real(nn(k)) c 425 continue c 430 continue 440 continue c c *****Same lasts if(ii .lt. 3)goto 471 do 470, k=3,ii do 460, l=1,(k-1) c do 455, t=1,(k-1) if(l .eq. t)goto 455 do 450, m=1,nn(k) cov12(l,k,t,k)=cov12(l,k,t,k)+(vv(1,k,l,m)- + area(1,l,k))*(vv(2,k,t,m)-area(2,t,k)) 450 continue c cov12(l,k,t,k)=(cov12(l,k,t,k)/real(nn(k)-1))/ + real(nn(k)) c 455 continue c 460 continue 470 continue c 471 continue c c *****Same middles do 490, k=2,(ii-1) do 485, l=1,(ii-2) if(l .ge. k)goto 485 c do 480, t=(k+1),ii do 475, m=1,nn(k) cov12(l,k,k,t)=cov12(l,k,k,t)+(vv(1,k,l,m)- + area(1,l,k))*(vv(2,k,t,m)-area(2,k,t)) 475 continue c cov12(l,k,k,t)=(cov12(l,k,k,t)/real(nn(k)-1))/ + real(nn(k)) c 480 continue c 485 continue 490 continue c c *****First of modality 1 equals last of modality 2 do 590, k=2,(ii-1) do 585, t=1,(ii-2) if(t .ge. k)goto 585 c do 580, l=(k+1),ii do 575, m=1,nn(k) cov12(k,l,t,k)=cov12(k,l,t,k)+(vv(1,k,l,m)- + area(1,l,k))*(vv(2,k,t,m)-area(2,k,t)) 575 continue c cov12(k,l,t,k)=(cov12(k,l,t,k)/real(nn(k)-1))/ + real(nn(k)) c 580 continue c 585 continue 590 continue c write(2,600) 600 format(9x,'Modality 1',12x,'Modality 2') write(2,601) 601 format(8x,'Area', 3x,'Var',6x,'Area', 3x,'Var',4x,'Cov') do 625, o=1,(ii-1) do 624, s=(o+1),ii write(2,620)o,s,area(1,o,s),var(1,o,s),area(2,o,s), + var(2,o,s),cov12(o,s,o,s) 620 format(i2,' v',i2,2x,f6.4,1x,f9.7,4x,f6.4,1x,f9.7,4x,f9.7) 624 continue 625 continue c c c Sum up the covariances c cov12s1=0.0 cov12s2=0.0 cov12s3=0.0 c do 640, o=1,ii do 639, s=1,ii do 638, t=1,ii do 637, u=1,ii c cov12s1=cov12s1+(real(nn(o)*nn(s))/ + denom1)*(real(nn(t)*nn(u))/denom1)*cov12(o,s,t,u) c cov12s2=cov12s2+(prev(s)*prev(o)/denom2) + *(prev(t)*prev(u)/denom2)*cov12(o,s,t,u) c cov12s3=cov12s3+(prev(o)*prev(s)/denom2) + *(prev(t)*prev(u)/denom2)*loss(o,s)*loss(t,u)* + cov12(o,s,t,u) c 637 continue 638 continue 639 continue 640 continue c write(2,641) 641 format(//) c c do 650, modal=1,2 ses1(modal)=sqrt(vars1(modal)+covars1(modal)) ses2(modal)=sqrt(vars2(modal)+covars2(modal)) ses3(modal)=sqrt(vars3(modal)+covars3(modal)) 650 continue c write(2,651) 651 format(20x,'Modality 1',13x,'Modality 2') write(2,652) 652 format(16x,'Summary', 2x,'SE',9x,'Summary', 2x,'SE',9x,'Cov') c write(2,655)summary1(1),ses1(1),summary1(2),ses1(2),(cov12s1) write(2,656)summary2(1),ses2(1),summary2(2),ses2(2),(cov12s2) write(2,657)summary3(1),ses3(1),summary3(2),ses3(2),(cov12s3) c c 655 format('Sample Prev ',2(f6.4,5x,f9.7,4x),f9.7) 656 format('Pop Prev ',2(f6.4,5x,f9.7,4x),f9.7) 657 format('Pop Prev w Loss ',2(f6.4,5x,f9.7,4x),f9.7) c 1000 continue stop end