mcn99a.f（Fortran）

＜mcn99a.f（Fortran）プログラムソース＞

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

　
＜プログラム処理結果サンプル＞