PROGRAM MCN99A
C     Fortran code for non-inferiority test, confidence interval 
C     for the difference in proportions for the paired-sample design
C
C                  August 1999 (version 2.00)
C
C                          Programmed by
C 
C                          Toshiro Tango
C                  Division of Theoretical Epidemiology
C                       Department of Epidemiology
C                     The Institute of Public Health
C              4-6-1 Shirokanedai, Minatoku, Tokyo, 108, Japan
C                        E-mail: tango@iph.go.jp
C
C ----------------------------------------------------------------
C                    (  README  FIRST )
C  1. Please read my paper carefully before using this program.
C  2. You are free to use this program for noncommercial purposes
C     even to modify it but please cite this original program
C 
C  Reference: 
C  Tango, T. Equivalence test and confidence interval for the difference 
C    in proportions for the paired-sample design. 
C                    Statistics in Medicine 17, 891-908 (1998).
C -----------------------------------------------------------------
C         
C      DIMENSION ORD(2),MRD(2),QHG(2)
	dimension Z(2)
      DIMENSION v1(100),v2(100),v3(100),v4(100),xx(2)
C      external random 
C      external ran_ini
C            
      real*4 EPS
      real*4 z
C	real*4 ord,mrd,qhg
C      REAL*4 RND
C      integer*4 init_val 
C      integer*4 random_seed
      integer*4 i, iii, stflag
      integer*4 v1,v2,v3,v4,a,b,c,d
      integer*4 npair,tra
      character*8 ti,tj
	external time
      common random_seed
C     input data file
	OPEN(7,FILE='d:\demo\mcn99a.dat')
C     output data.file
 	OPEN(8,FILE='d:\demo\mcn99a.out')           
C     COMMON /RANDOM/ RANDOM_SEED      
C
      call time( ti )
      write(8,*) ti
      read(7,*)npcnt,zval,delta,npair
C
C     npcnt: percentage, for example, = 95 for 95 % CI
C     zval : normal deviate for npcnt, for example, =1.960 for 95%
C     delta: delta value for equivalence test, for example, = 0.1
C     npair: number of pairs of (a,b,c,d) for calculation
C     see input data file mcn99a.dat
C
      do 100 i=1,npair
      read(7,*) v1(i),v2(i),v3(i),v4(i)
C     these four v's are (a,b,c,d)
  100 continue
C
      DO 999 III=1,NPAIR
C
      a=v1(iii)
      b=v2(iii)
      c=v3(iii)
      d=v4(iii)
      n=a+b+c+d
	xn=n
C      write(6,246)a,b,c,d,a+b+c+d
C  246 format(1H ,5F5.0)
C
      EPS=0.00001
      z(1)= -zval
      z(2)= -z(1)
C
C     Test for equivalence
C
      S=delta
      if((s .lt. 0.001) .and. (b .eq. 0) .and. (c .eq. 0))then
            t4=0
            goto 430
      endif
  430 tra=0
      CALL MCNE(S,float(N),float(B),float(C),T4)
      CALL TAN(ZVAL,A,B,C,D,XX,TRA,STFLAG)
C
      write(8,833)a,b,c,d,n,xx(1),xx(2),
     *STFLAG,T4,tra
C
C     xx(1): lower bound of CI
C     xx(2): upper bound of CI
C     STFLAG: =0, for convergence; =1, for not converged !
C     T4: the proposed test statistics for equivalence test
C     tra=0 if b >= c; =1 for if b < c,   just for check
C     see example output file: mcn99a.out
C 
  833 format(1H ,5I4,2F10.5,I5,F10.5,I5)
  999 continue
C
      call time( tj )
      write(8,*) tj
      STOP 
      END
C
      SUBROUTINE TAN(ZVAL,A,B,C,D,XX,TRA,STFLAG)
      DIMENSION XX(2),SC(2),Z(2)
      INTEGER*4 A,B,C,D,N,T,tra,stflag,flag
      N=A+B+C+D
      XN=N
      EPS=0.00001
      z(1)= -zval
      z(2)= -z(1)
      tra=0
      if(b .lt. c)then
         t=c
         c=b
         b=t
         tra=1
      endif
      ns=1
      ne=2
      if((b .eq. 0) .and. (c .eq. 0))then
         ddd=z(1)**2/(xn+z(1)**2)
         sc(1)=ddd
         sc(2)=-ddd
         goto 682
      endif
      if((a .eq. 0) .and. (c .eq. 0) .and.
     * (d .eq. 0))then
         ns=2
         ne=2
      endif
C   
      do 600 j=ns,ne
      if((b .gt. 0) .and. (c .gt. 0))then
         flag=1
         goto 520
      endif
      flag=2
  520 if(flag .eq. 1)then
          bh=b
	  ch=c
          x0=(bh-ch)/xn-z(j)/xn*sqrt( ch+bh-(ch-bh)**2/xn )
	endif
      if((abs(x0) .gt. 1) .or. (flag .eq. 2))x0=0.998
      ZZ=Z(J)
      call comp(ZZ,EPS,X0,X1,float(B),float(C),float(N),SS,stflag)
      SC(J)=SS
  600 continue
C
      if((a .eq. 0) .and. (c .eq. 0) .and.
     * (d .eq. 0))sc(1)=1
  682 xx(1)=sc(2)
      xx(2)=sc(1)
      if(tra .eq. 1)then
	  xx(1)=-sc(1)
	  xx(2)=-sc(2)
      endif
      if(tra .eq. 1)then
        t=c
        c=b
        b=t
      endif
      RETURN
      END

C    
      SUBROUTINE COMP(ZZ,EPS,X0,X1,GB,GC,N,SS,stflag)
      REAL*4 ZZ,X0,X1,GB,GC,N,SS,Y0,Y1,W,V,X,EPS,T3,S,H0,H1
	integer*4 stflag
      X1=X0*0.95
      ncount=0
      stflag=0
  300 S= -X0
      ncount=ncount+1
      if(ncount .gt. 2000)then
         stflag=1
         goto 310
      endif
      CALL MCNE(S,N,GB,GC,T3)
      Y0=T3
      S= -X1
      CALL MCNE(S,N,GB,GC,T3)
      Y1=T3
      H0=ABS(Y0+ZZ)
      H1=ABS(Y1+ZZ)
      IF(ABS((Y0-Y1)/Y0) .LE. EPS) GOTO 310
      IF( H0 .GT. H1) THEN 
        W=X0
        X0=X1
        X1=W
        V=Y0
        Y0=Y1
        Y1=V
      ENDIF
      X=X0+(X1-X0)*(ZZ-Y0)/(Y1-Y0)
      IF(ABS((X-X0)/X0) .LT. EPS) GOTO 310
      X0=X1
      X1=X
C       WRITE(6,*)X0,X1,Y0,Y1
      GOTO 300
  310 SS=X
      RETURN
      END      
C
      SUBROUTINE MCNE(S,N,GB,GC,T3)
      REAL*4 VA,VB,VC,VQ21,VQ12,T3,N,S,GB,GC
      VA=2*N
      VB= -(2*N-GB+GC)*S -GB-GC
      VC= GC*S*(1.0+S)
      VQ21=(SQRT(VB*VB-4*VA*VC)-VB)/2.0/VA
      VQ12=VQ21-S
      T3=(GB-GC+N*S)/SQRT( N*(VQ21+VQ12-S*S) )
      RETURN
      END