c	Interactive FORTRAN program, DESIGNROC, to compute sample size/power 
c	using formulae in Obuchowski and McClish, Stat in Med 1997; 16: 
c	1529-1542. This program is for comparing two diagnostic tests, 
c	paired or unpaired design, using ROC area, partial area, or Sen 
c	at fixed FPR.   
c	
c	written by Nancy A. Obuchowski, The Cleveland Clinic Foundation
c	Cleveland, Ohio.
c
	implicit real(a-z)
	double precision parta,cum,pd
	integer index,pwss,new
	character*3 answ
c
	pi=3.1415927
	ff=0.0
	gg=0.0
	new=0
c
	print *,'Make your selection of desired index:'
	print *,'area under full curve [1]' 
	print *,'partial area under curve [2]'
	print *,'sensitivity at fixed FPR [3]'
	read *,index
c
c
c
	print *,'Enter the binormal parameters A and B for the first test'
	read *,aa1,bb1
c
	print *,'Enter the desired ratio of normal/abnormal subjects'
	read *,rr
c
c
	expr1=exp((-aa1*aa1)/(2.0*(1.0+bb1*bb1)))
	expr2=1.0+bb1*bb1
c
c
	if(index .eq. 1)then
		gg1=-expr1*aa1*bb1/sqrt(2.0*pi*expr2**3)
		ff1=expr1/sqrt(2.0*pi*expr2)
	endif
c
c
 10	continue
c	**************************************************
	if(index .eq. 2)then
	print *,'enter the FPR range of interest'
	read *,FPR1,FPR2
	c1=zscore(FPR1)
	c2=zscore(FPR2)
	if(FPR1 .lt. 0.5)c1=-c1
	if(FPR2 .lt. 0.5)c2=-c2
	if(FPR1 .le. 0.0)c1=-6.0
	if(FPR2 .ge. 1.0)c2=6.0
c	print *,'c1,c2',c1,c2
	c1p=(c1+aa1*bb1/expr2)*sqrt(expr2)
	c2p=(c2+aa1*bb1/expr2)*sqrt(expr2)
	c1pp=(c1p**2)/2.0
	c2pp=(c2p**2)/2.0
c	print *,'c1p,c2p,c1pp,c2pp',c1p,c2p,c1pp,c2pp
	expr3=(1.0-pnorm(c2p))-(1.0-pnorm(c1p))
	expr4=exp(-c1pp)-exp(-c2pp)
	ff1=expr1/sqrt(2.0*pi*expr2)*expr3
	gg1=expr1/(2.0*pi*expr2)*expr4-aa1*bb1*expr1/
     + sqrt(2.0*pi*expr2**3)*expr3
	endif
c
c	****************************************************
c
	if(index .eq. 3)then
	print *,'enter fixed FPR of interest'
	read *,FPR
	ff1=1.0
	gg1=zscore(1.0-FPR)
	endif
c
c	print *,'f1,g1: ',ff1,gg1
c
c
c
	if(new .eq. 1)goto 12
c
	print *,'Enter the binormal parameters A and B for the second test'
	read *,aa2,bb2
c
c
	expr11=exp((-aa2*aa2)/(2.0*(1.0+bb2*bb2)))
	expr12=1.0+bb2*bb2
c
c
	if(index .eq. 1)then
		gg2=-expr11*aa2*bb2/sqrt(2.0*pi*expr12**3)
		ff2=expr11/sqrt(2.0*pi*expr12)
	endif
c
 12	continue
c
c	**************************************************
	if(index .eq. 2)then
	c1p=(c1+aa2*bb2/expr12)*sqrt(expr12)
	c2p=(c2+aa2*bb2/expr12)*sqrt(expr12)
	c1pp=(c1p**2)/2.0
	c2pp=(c2p**2)/2.0
	expr3=(1.0-pnorm(c2p))-(1.0-pnorm(c1p))
	expr4=exp(-c1pp)-exp(-c2pp)
	ff2=expr11/sqrt(2.0*pi*expr12)*expr3
	gg2=expr11/(2.0*pi*expr12)*expr4-aa2*bb2*expr11/
     + sqrt(2.0*pi*expr12**3)*expr3
	endif
c
c	****************************************************
c
	if(index .eq. 3)then
	ff2=ff1
	gg2=gg1
	endif
c
c	print *,'f2,g2: ',ff2,gg2
c
c
	if(new .eq. 1)goto 15
	print *,'What is A and B under the null hypothesis?'
	read *,aao,bbo
c
	expr111=exp((-aao*aao)/(2.0*(1.0+bbo*bbo)))
	expr112=1.0+bbo*bbo
c
c
	if(index .eq. 1)then
		ggo=-expr111*aao*bbo/sqrt(2.0*pi*expr112**3)
		ffo=expr111/sqrt(2.0*pi*expr112)
	endif
c
 15	continue
c
c	**************************************************
	if(index .eq. 2)then
	c1p=(c1+aao*bbo/expr112)*sqrt(expr112)
	c2p=(c2+aao*bbo/expr112)*sqrt(expr112)
	c1pp=(c1p**2)/2.0
	c2pp=(c2p**2)/2.0
	expr113=(1.0-pnorm(c2p))-(1.0-pnorm(c1p))
	expr114=exp(-c1pp)-exp(-c2pp)
	ffo=expr111/sqrt(2.0*pi*expr112)*expr113
	ggo=expr111/(2.0*pi*expr112)*expr114-aao*bbo*expr111/
     + sqrt(2.0*pi*expr112**3)*expr113
	endif
c
c	****************************************************
c
	if(index .eq. 3)then
	ffo=ff1
	ggo=gg1
	endif
c
c
c
	if(new .eq. 1)goto 20
	print *,'Enter the correlations in the underlying distribution for'
	print *,'the test results of nondiseased and diseased subjects'
	read *,rn,ra
	new=1
c
 20	continue
c
c
c	Variance under the null hypothesis of no difference is:
c
	varo=ffo*ffo*(1.0+bbo*bbo/rr+aao*aao/2.0)+ggo*ggo*
     + (bbo*bbo*(1.0+rr))/(2.0*rr)
c
	varob=ffo*ffo*(1.0+bbo*bbo/rr+aao*aao/2.0)+ggo*ggo*
     + (bbo*bbo*(1.0+rr))/(2.0*rr)+2.0*ffo*ggo*(aao*bbo)/2.0
c
	covo=ffo*ffo*(ra+(rn*bbo*bbo)/rr+(ra*ra*aao*aao)/2.0)+
     + ggo*ggo*((rn*rn*bbo*bbo)/(2.0*rr)+(ra*ra*bbo*bbo)/2.0)+
     + 2.0*ffo*ggo*(ra*ra*aao*bbo)/2.0
c	print *,'varo,covo',varo,covo
c
	varnull=2.0*varo-2.0*covo
c	print *,'variance of index under null ',varnull
	varnullb=2.0*varob-2.0*covo
c	print *,'variance of index under direct estimation', varnullb
c
c
c	Variance under the alternative hypothesis is:
c
	vara1=ff1*ff1*(1.0+bb1*bb1/rr+aa1*aa1/2.0)+gg1*gg1*
     + (bb1*bb1*(1.0+rr))/(2.0*rr)
c
	vara1b=ff1*ff1*(1.0+bb1*bb1/rr+aa1*aa1/2.0)+gg1*gg1*
     + (bb1*bb1*(1.0+rr))/(2.0*rr)+2.0*ff1*gg1*(aa1*bb1)/2.0
c
	vara2=ff2*ff2*(1.0+bb2*bb2/rr+aa2*aa2/2.0)+gg2*gg2*
     + (bb2*bb2*(1.0+rr))/(2.0*rr)
c
	vara2b=ff2*ff2*(1.0+bb2*bb2/rr+aa2*aa2/2.0)+gg2*gg2*
     + (bb2*bb2*(1.0+rr))/(2.0*rr)+2.0*ff2*gg2*(aa2*bb2)/2.0
c
	cova=ff1*ff2*(ra+(rn*bb1*bb2)/rr+(ra*ra*aa1*aa2)/2.0)+
     + gg1*gg2*((rn*rn*bb1*bb2)/(2.0*rr)+(ra*ra*bb1*bb2)/2.0)+
     + ff1*gg2*(ra*ra*aa1*bb2)/2.0+ff2*gg1*(ra*ra*aa2*bb1)/2.0
c
	varalt=vara1+vara2-2.0*cova
c	print *,'variance of index under alternative ',varalt
	varaltb=vara1b+vara2b-2.0*cova
c
c
c	Compute the difference to be detected
c
	if(index .eq. 1)then
c	print *,' ... working'
	ac1=-6.0
	ac2=6.0
	aa=aa1
	bb=bb1
	parta=0.0
	DO 25,ni=ac1,ac2,0.005
	cum=1.0-PNORM(aa+bb*ni)
	pd=(1.0/SQRT(2.0*pi))*(EXP(-0.5*(ni*ni)))
	parta=parta+(cum*pd*0.005)
 25	CONTINUE
	area1=parta
c
	aa=aa2
	bb=bb2
	parta=0.0
	DO 26,ni=ac1,ac2,0.005
	cum=1.0-PNORM(aa+bb*ni)
	pd=(1.0/SQRT(2.0*pi))*(EXP(-0.5*(ni*ni)))
	parta=parta+(cum*pd*0.005)
 26	CONTINUE
	area2=parta
	diff=abs(area1-area2)
c	print *,diff
	endif
c
c
	if(index .eq. 2)then
c	print *,' ... working'
	aa=aa1
	bb=bb1
	parta=0.0
	DO 30,ni=c1,c2,0.0005
	cum=1.0-PNORM(aa+bb*ni)
	pd=(1.0/SQRT(2.0*pi))*(EXP(-0.5*(ni*ni)))
	parta=parta+(cum*pd*0.0005)
 30	CONTINUE
	area1=parta
c
	aa=aa2
	bb=bb2
	parta=0.0
	DO 31,ni=c1,c2,0.0005
	cum=1.0-PNORM(aa+bb*ni)
	pd=(1.0/SQRT(2.0*pi))*(EXP(-0.5*(ni*ni)))
	parta=parta+(cum*pd*0.0005)
 31	CONTINUE
	area2=parta
	diff=abs(area1-area2)
c	print *,diff
	endif
c
c
	if(index .eq. 3)then
	area1=aa1+bb1*gg1
	area2=aa2+bb2*gg2
	diff=abs(area1-area2)
	endif	
c
c
	print *,'For one tailed test, enter the type one error rate;'
	print *,'for two tailed test, enter area of one of the tails'
	read *,alpha
c
 100	continue
	print *,'Do you want to estimate power [1] or sample size [2]?'
	read *,pwss
c
	if(pwss .eq. 2)then
	print *,'Enter the desired power'
	read *,power
c
c
	zalpha=zscore(alpha)
	zbeta=zscore(1.0-power)
c
c
c
	na=((zalpha*sqrt(varnull)+zbeta*sqrt(varalt))**2)/
     + (diff*diff)
	nn=na*rr
	print *,'*********************************************************'
	print *,'*********************************************************'
	print *,'THE REQUIRED SAMPLE SIZE IS:'
	print *,na,' abnormals and ',nn,' normals'
c
	print *,'*********************************************************'
	print *,'In the SPECIAL CASE where the ROC curve will be estimated'
	print *,'directly from the test results, the following sample'
	print *,'size approximation may be reasonable'
	na=((zalpha*sqrt(varnullb)+zbeta*sqrt(varaltb))**2)/
     + (diff*diff)
	nn=na*rr
c	print *,'required sample size is:'
	print *,na,' abnormals and ',nn,' normals'
	print *,'********************************************************'
	endif
c	
c
 	if(pwss .eq. 1)then
	zalpha=zscore(alpha)
	print *,'Enter the number of diseased subjects'
	read *,na
	zbeta=(sqrt(na)*diff-zalpha*sqrt(varnull))/sqrt(varalt)
	power=pnorm(zbeta)
	print *,'*********************************************************'
	print *,'*********************************************************'
	print *,'THE ESTIMATED POWER IS ',(1.0-power)
c
	print *,'*********************************************************'
	print *,'In the SPECIAL CASE where the ROC curve will be estimated'
	print *,'directly from the test results, the following power'
	print *,'estimate may be reasonable'
	zbeta=(sqrt(na)*diff-zalpha*sqrt(varnullb))/sqrt(varaltb)
	power=pnorm(zbeta)
	print *,(1.0-power)
	print *,'********************************************************'

	endif
c
	print *,'Would you like to keep everything constant but change' 
	print *,'the sample size or power?'
	read *,answ
	if(answ .eq. 'YES' .or. answ .eq. 'yes' .or. answ .eq. 'y'
     + .or. answ .eq. 'Y')goto 100
c
	if(index .eq. 1)goto 110
c
	print *,'Would you like to keep everything constant but change'
	print *,'the FPR of interest?'
	read *,answ
	if(answ .eq. 'YES' .or. answ .eq. 'yes'.or. answ .eq. 'y'
     + .or. answ .eq. 'Y')goto 10
c
 110	continue
	end
c
c	********************************************
	function zscore(pvalue)
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
	if(pvalue .lt. 1.0)goto 30
	ifault=1
	return
 30	continue
	if(pvalue .gt. 0)goto 40
	ifault=2
	return
 40	continue
	if(pvalue .le. 0.50)goto 50
	sign=-1.0
	pvalue=1.0-pvalue
 50	continue
	q=pvalue
	t=sqrt(alog(1/Q**2))
	den=1+d1*t+d2*(t**2)+d3*(t**3)
	num=c0+c1*t+c2*(t**2)
	zscore=(t-num/den)*sign
	return
	end
c
c
c	************************************************
	function pnorm (zscore)
c	************************************************
	implicit real (a-z)
	r=0.2316419
	pi=3.14159268
	b1=0.31938153
	b2=-0.356563782
	b3=1.781477937
	b4=-1.821255978
	b5=1.330274429
	sign=1.0
	ref=0.0
	if(zscore .ge. 0.0)goto 50
	sign=-1.0
	ref=1.0
	zscore=-1.0*zscore
 50	continue
	x=zscore
	t=1/(1+r*x)
	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))
	pnorm=sign*(func*poly)+ref
	return
	end