Modify fortran tests for debugging chi-square issue

Author: viking

misc/fortran-tests/Makefile

@@ -1,16 +1,7 @@
 FC = gfortran -std=legacy
 
-ps-test: zcrvalue.o pone.o ssize.o main.o
-	$(FC) -o ps-test zcrvalue.o pone.o ssize.o main.o
+ps-test: zcrvalue.o pone.o ssize.o iprelris.o chisqsiz.o bisec.o fishersi.o main.o
+	$(FC) -o $@ $^
 
-zcrvalue.o: zcrvalue.for
-	$(FC) -c -o zcrvalue.o zcrvalue.for
-
-pone.o: pone.for
-	$(FC) -c -o pone.o pone.for
-
-ssize.o: ssize.for
-	$(FC) -c -o ssize.o ssize.for
-
-main.o: main.for
-	$(FC) -c -o main.o main.for
+%.o: %.for
+	$(FC) -c -o $@ $<

misc/fortran-tests/bisec.for

@@ -0,0 +1,114 @@
+	SUBROUTINE BISEC(Y1_ARG,Y2_ARG,E_ARG,E1_ARG,F,X,ERROR)
+C
+C-Description-----------------------------------------------------------
+	IMPLICIT NONE
+C
+C  Function:
+C	bisection
+C	This procedure evaluates a function at the end-points of a real
+C	interval.  An error exit is taken if there is no change of sign.
+C	Finds a root by iterated bisection and evaluation at the midpoint,
+C	halting if either the value of the function is less than E_ARG or
+C	two successive approximations of the root differ by less than
+C	E1_ARG.
+C  Usage:
+C	
+C  Arguments:
+C   i	Y1_ARG	Lower end of interval
+C   i	Y2_ARG  Upper end of interval
+C   i	E_AGR	If value of function <= E_ARG then done
+C   i	E1_ARG  If two successive approximations differ <= E1_ARG then done
+C   i	F	The function to evaluate
+C   o	X	Argument to function that produces root if one is found
+C   o	ERROR	0 if solution found
+C		1 if no solution is found
+C  Notes:
+C   .	Reference:
+C	"Bisection Routine", _Collected Algorithms from ACM_, vol I,
+C	Algorithm 4, 1980.
+C
+C-Declarations----------------------------------------------------------
+C
+C
+C  Arguments
+C
+	REAL Y1_ARG,Y2_ARG,E_ARG,E1_ARG,X
+	INTEGER ERROR
+C
+C  Functions
+C
+	REAL F
+C
+C  Locals
+C
+        REAL F1,FX
+	REAL Y1,Y2,E,E1
+	INTEGER I,J,K
+C
+C-Code------------------------------------------------------------------
+C
+C    Assume success.
+C
+	ERROR=0
+C
+C    Copy arguments to local variables.
+C
+	Y1=Y1_ARG
+	Y2=Y2_ARG
+	E=E_ARG
+	E1=E1_ARG
+C
+C    Initialize values.
+C
+10	CONTINUE
+	I=1
+	J=1
+	K=1
+	X=Y2
+C
+C    Evaluate the function at X.  If result is <= E then done.
+C
+20	FX=F(X)
+	IF (ABS(FX).LE.E) RETURN
+	GO TO (30,40) I
+30	CONTINUE
+	I=2
+	F1=FX
+	X=Y1
+	GO TO 20
+40	CONTINUE
+	IF (FX*F1.GT.0) GO TO (1090,70) J
+	GO TO (50,90) K
+50	CONTINUE
+	J=2
+	K=2
+C
+C    Calculate a new midpoint.
+C
+60	CONTINUE
+	X=(Y1+Y2)/2.
+	GO TO 20
+C
+C    Set interval end to be current midpoint.
+C
+70	CONTINUE
+	Y2=X
+C
+C    If two successive iteration differ by <= E1 then end.
+C
+80	CONTINUE
+	IF (ABS(Y1-Y2).GE.E1) GO TO 60
+	RETURN
+C
+C    Set interval end to be current midpoint.
+C
+90	CONTINUE
+	Y1=X
+	GO TO 80
+1090	CONTINUE
+C
+C    Set an error flag if unable to find a root.
+C
+	ERROR=1
+	RETURN
+	END

misc/fortran-tests/chisqsiz.for

@@ -0,0 +1,42 @@
+	REAL FUNCTION CHISQSIZE(P1)
+	IMPLICIT NONE
+C
+C    This function returns the sample size associated with P1 and the
+C    other input parameters.  This formula should be used for studies that
+C    are analyzed with an uncorrected chi-squared test.
+C
+C    Arguments
+C
+	REAL P1
+C
+C    Functions
+C
+	REAL ZCRVALUE
+C
+C    Locals
+C
+	REAL ZALPHA,ZBETA,P,Q
+C
+C    Common from main
+C
+	REAL ALPHA,BETA,P0,N,M,Q0,Q1
+	COMMON/IPCOM/ALPHA,BETA,P0,N,M,Q0,Q1
+C
+C Code:
+	ZALPHA=ZCRVALUE(ALPHA/2.)
+	ZBETA=ZCRVALUE(BETA)
+	P=(N*P1+M*N*P0)/(M*N+N)
+        Q=1.-P
+	Q1=1.-P1
+	Q0=1.-P0
+	CHISQSIZE=
+     .	(
+     .	(ZALPHA*SQRT((1.+1./M)*P*Q) + ZBETA*SQRT((P0*Q0)/M+(P1*Q1)))**2
+     .     /
+     .  (P0-P1)**2
+     .	) - N
+C
+C    Done.
+C
+	RETURN
+	END

misc/fortran-tests/fishersi.for

@@ -0,0 +1,42 @@
+	REAL FUNCTION FISHERSIZ(P1)
+	IMPLICIT NONE
+C
+C    This function returns the sample size associated with P1 and the
+C    other input parameters.  This formula should be used for studies that
+C    are analyzed with Fisher's Exact test.
+C
+C    Arguments
+C
+	REAL P1
+C
+C    Functions
+C
+	REAL ZCRVALUE
+C
+C    Locals
+C
+	REAL ZALPHA,ZBETA,P,Q,NPRIME
+C
+C    Common from main
+C
+	REAL ALPHA,BETA,P0,N,M,Q0,Q1
+	COMMON/IPCOM/ALPHA,BETA,P0,N,M,Q0,Q1
+C
+C Code:
+	ZALPHA=ZCRVALUE(ALPHA/2.)
+	ZBETA=ZCRVALUE(BETA)
+	P=(N*P1+M*N*P0)/(M*N+N)
+        Q=1.-P
+	Q1=1.-P1
+	Q0=1.-P0
+	NPRIME=
+     . (ZALPHA*SQRT((1.+1./M)*P*Q) + ZBETA*SQRT((P0*Q0)/M+(P1*Q1)))**2
+     .     /
+     .  (P0-P1)**2
+	FISHERSIZ=
+     . (NPRIME*(1.+SQRT(1.+2.*(M+1.)/(NPRIME*M*ABS(P0-P1))))**2/4.) - N
+C
+C    Done.
+C
+	RETURN
+	END

misc/fortran-tests/iprelris.for

@@ -0,0 +1,144 @@
+
+	SUBROUTINE IPRELRISK
+     . (XALPHA,POWER,XP0,XN,XM,UORF,P1L,ODDSL,RL,P1H,ODDSH,RH)
+!MS$ATTRIBUTES REFERENCE :: P1L
+!MS$ATTRIBUTES REFERENCE :: ODDSL
+!MS$ATTRIBUTES REFERENCE :: RL
+!MS$ATTRIBUTES REFERENCE :: P1H
+!MS$ATTRIBUTES REFERENCE :: ODDSH
+!MS$ATTRIBUTES REFERENCE :: RH
+	IMPLICIT NONE
+
+	REAL*4 XALPHA,POWER,XP0,XN,XM,P1L,ODDSL,RL,P1H,ODDSH,RH
+	INTEGER*4 UORF
+C
+C	Relative risk CALCULATIONS FOR INDEPENDENT PROPORTIONS
+C	DERIVE Relative risk FOR EACH GROUP IN A TRIAL OF TWO TREATMENTS 
+C	OR A CASE-CONTROL STUDY.
+C	NORMAL APPROXIMATION ASSUMED
+C
+C       Power formula from:  Schlesselman:  Case-control Studies: Design,
+C       Conduct, Analysis.  New York: Oxford U. Press 1982: 144-152.
+C		See also:  FRIEDMAN, FURBERG AND DEMETS: 
+C		FUNDAMENTALS OF CLINICAL TRIALS. BOSTON: WRIGHT PSG. 1982:75.
+C
+C INPUT:
+C	ALPHA	Type I error probability (two tailed)
+C	POWER	The desired statistical power
+C	P0	True probability of exposure in control group when CASECTL=Y
+C		True probability of event in control group when CASECTL=N
+C	N	Number of case patients
+C	M	Ratio of control to case
+C	CASECTRL 1 for case-control studies
+C			2 for prospective studies
+C	UORF =1 for uncorrected chi-square test
+C        =2 for Fisher's exact test
+C OUTPUT:
+C	RELRISK	Relative risk of cases relative to controls that can be
+C		detected with power POWER (when casectl=N)
+C	ODDSRATIO Oddsratio of cases to controls that can be detected
+C		with power POWER (when casectl=Y)
+C	P1	True probability of exposure in case group that can be
+C		detected with power POWER when CASECTL=Y
+C		True probability of event in experimantal group that can
+C		be detected with power POWER when CASECTL=N
+C
+C
+C  Globals
+C
+C    Declare variables in common to be made available to the function.
+C
+	REAL ALPHA,BETA,P0,N,M,Q0,Q1
+	COMMON/IPCOM/ALPHA,BETA,P0,N,M,Q0,Q1
+C
+C    Declare functions.
+C
+	EXTERNAL CHISQSIZE
+	REAL CHISQSIZE
+	EXTERNAL FISHERSIZ
+	REAL FISHERSIZ
+C
+C    Locals
+C
+	INTEGER ERRORL,ERRORH
+	REAL E1,EPS,Y1,Y2,Q1L,Q1H
+C
+C Code:
+C
+	ALPHA=XALPHA
+	BETA=1.-POWER
+	P0=XP0
+	N=XN
+	M=XM
+C
+C    Set error tolerance values.
+C
+	EPS=.0001*MIN(P0,1.-P0)
+	E1=.0001*MIN(P0,1.-P0)
+C
+C    Use bisection with appropriate routine for Chi-squared or Fisher's
+C    exact test.
+C
+	IF (UORF.EQ.1) THEN
+C
+C    Set initial end points for lower end point.
+C
+	    Y1=0.0
+	    Y2=P0-EPS
+C
+C    Solve the equation for the lower solution.
+C
+	    CALL BISEC(Y1,Y2,EPS,E1,CHISQSIZE,P1L,ERRORL)
+C
+C    Set initial end points for high end point.
+C
+	    Y1=P0+EPS
+	    Y2=1.
+C
+C    Solve the equation for the upper solution.
+C
+	    CALL BISEC(Y1,Y2,EPS,E1,CHISQSIZE,P1H,ERRORH)
+	ELSE
+C
+C    Case of Fisher's exact test.
+C
+	    Y1=0.0
+	    Y2=0.9999*P0
+	    CALL BISEC(Y1,Y2,EPS,E1,FISHERSIZ,P1L,ERRORL)
+	    Y1=0.9999*P0+0.0001
+	    Y2=1.
+	    CALL BISEC(Y1,Y2,EPS,E1,FISHERSIZ,P1H,ERRORH)
+	END IF
+C
+C    Calculate the relative risk and odds ratio for low end.
+C    21-Jan-1992 - Added calculation of Q1L to fix calculation of ODDSL.
+C
+	IF (ERRORL.EQ.0) THEN
+		RL=P1L/P0
+		Q1L=1.-P1L
+		ODDSL=(P1L*Q0)/(P0*Q1L)
+	ELSE
+		P1L=0
+		RL=0
+		Q1L=0
+		ODDSL=0
+	END IF
+C
+C    Calculate the relative risk and odds ratio for high end.
+C    21-Jan-1992 - Added calculation of Q1H to fix calculation of ODDSH.
+C
+	IF (ERRORH.EQ.0) THEN
+		RH=P1H/P0
+		Q1H=1.-P1H
+		ODDSH=(P1H*Q0)/(P0*Q1H)
+	ELSE
+		P1H=0
+		RH=0
+		Q1H=0
+		ODDSH=0
+	END IF
+C
+C	Done.
+C
+	RETURN
+	END

misc/fortran-tests/main.for

@@ -1,16 +1,17 @@
 C vim: set noexpandtab:
 	PROGRAM MAIN
 	IMPLICIT NONE
-	REAL ALPHA, BETA, PHI, P0, PSI, XN
-	INTEGER M, ERR
+	REAL*4 ALPHA,POWER,P0,N,M,P1L,ODDSL,RL,P1H,ODDSH,RH
+	INTEGER*4 UORF
 
 	ALPHA = 5D-2
-	BETA = 2D-1
-	PHI = 7D-1
-	P0 = 3D-1
-	M = 1
-	PSI = 2.
-	CALL SSIZE(ALPHA, BETA, PHI, P0, M, PSI, XN, ERR)
+	POWER = 8D-1
+	P0 = 1D-1
+	N = 325
+	M = 9
+	UORF = 1
 
-	PRINT *, "SSIZE(0.05, 1-0.8, 0.7, 0.3, 1, 2) =", XN, ERR
+	CALL IPRELRISK(ALPHA,POWER,P0,N,M,UORF,P1L,ODDSL,RL,P1H,ODDSH,RH)
+
+	PRINT *, "IPRELRISK(0.05, 0.8, 0.1, 325, 9, 1) =", P1L, P1H
 	END