Jung SH, Chow SC. On sample size calculation for comparing survival curves under general hypotheses testing, J of Biopharmaceutical Statistics, 2012; 22: 485-495.
The standard log-rank test is commonly used to test if an experimental arm has a longer survival distribution than a control arm, while the non-inferiority log-rank test is to test if an experimental arm has a survival distribution no better than but similar to that of a control arm. In this paper, the authors extend the use of the log-rank test for clinical trials with various types of non-conventional study objectives and propose its sample size calculation under general null and alternative hypotheses. The Fortran program calculate the sample size of randomized clinical trials comparing the survival distributions under various types of null and alternative hypotheses.
Fortran Begins ===================>
c sampLe size for noninferiority test using the Log-rank test
c Group 1=standard Rx; 2=experimentaL Rx
c Input: Lamda(1)=hazard rate for standard Rx, deLta=hz(2)/hz(1)
c r=accruaL rate, b=foLLow-up period, zaLpha, zbeta
c Output: a=accruaL period, n=sampLe size, events=required # events
integer n0(2)
common n(2),x(2,1000),d(2,1000),hz(2),p(2),
c deLta0,deLta1,b,r,za,zb,a0
idum=-1234
cccccccccccccccccccccccc
c Input parameters
ccccccccccccccccccccccccc
write(6,*)
write(6,*) '*** Input'
c # simulations to evaluate the estimated N.
c If you donot want simuLations, enter nsampLe=0.
nsim=10000
c st=survivaL probabiLity at time t0
t0=3
c Under H0
st1=0.85
hz(1)=-aLog(st1)/t0
delta0=3.84
hz20=delta0*hz(1)
write(6,*) t0,'-yr surv under H0=',exp(-hz(1)*t0),exp(-hz20*t0)
c st20=0.4
c write(6,*) t0,'-yr surv under H0=',st1,st20
c hz20=-aLog(st20)/t0
c deLta0=hz20/hz(1)
c Under H1
c st2=0.7
c delta1=2
c hz(2)=delta1*hz(1)
st2=0.85
hz(2)=-aLog(st2)/t0
write(6,*) t0,'-yr surv under H1=',exp(-hz(1)*t0),exp(-hz(2)*t0)
deLta1=hz(2)/hz(1)
write(6,101) hz(1),hz20,hz(1),hz(2)
101 format(2x,'Haz rate=',2f7.4,' under H0, ',2f7.4,' under H1')
write(6,102) deLta0,deLta1
102 format(2x,'HR=',f7.4,' under H0, ',f7.4,' under H1')
c HR: deLta0 (>1) under H0, deLta1 (<deLta0) under H1
c r=accruaL rate
r=60
c f=foLLow-up period
b=3
c p=aLLocation proportion for group 1
p(1)=.9
p(2)=1-p(1)
write(6,103) r,p(1),p(2)
103 format(2x,'Accrual rate=',f7.1,' alloc prop=',2f6.2)
c Type I and II errors
aLpha=0.05
beta=0.19
write(6,104) alpha,1-beta
104 format(2x,'Alpha=',f6.3,' 1-beta=',f6.3)
cccccccccccccccccccccccccccc
write(6,*)
write(6,*) '*** Output'
za=znorm(aLpha)
zb=znorm(beta)
c Bisection
a1=0.01
a2=10
f1=fun(a1)
f2=fun(a2)
if(f1*f2.Lt.0) go to 10
do k=1,10
a1=a2
a2=a1+10
f1=f2
f2=fun(a2)
if(f1*f2.Lt.0) go to 10
enddo
pause 'something wrong - 1'
10 continue
do k=1,20
a3=(a1+a2)/2.
f3=fun(a3)
if(r*abs(a1-a3).Lt.1.and.abs(f3).Lt.0.001) go to 12
if(f1*f3.Lt.0) then
a2=a3
f2=f3
eLse
a1=a3
f1=f3
endif
enddo
pause 'something wrong - 2'
12 continue
a0=a2
write(6,105) a0,b,a0+b
105 format(2x,'acc per=',f5.2,' f/u per=',f5.1,' (total=',f5.2,')')
nn=a0*r+1
n0(1)=nn*p(1)
n0(2)=nn*p(2)
nres=nn-n0(1)-n0(2)
d1=n0(1)*(1-exp(-hz(1)*b)*(1-exp(-hz(1)*a0))/a0/hz(1))
d2=(nn-n0(1))*(1-exp(-hz(2)*b)*(1-exp(-hz(2)*a0))/a0/hz(2))
dd=d1+d2
write(6,106) nn,n0(1),nn-n0(1)
106 format(2x,'N=',i5,' (n1 n2=',2i5,')')
write(6,107) dd,d1,d2
107 format(2x,'# events=',f7.1,' (nd1 nd2=',2f8.1,')')
c start simuLation
do 888 itype=2,1,-1
c To check empirical type I error, set hz(2)=hz20
if(itype.eq.1) hz(2)=hz20
pwr=0
do 999 isim=1,nsim
do k=1,2
n(k)=n0(k)
enddo
do l=1,nres
if(ran2(idum).Le.p(1)) then
n(1)=n(1)+1
eLse
n(2)=n(2)+1
endif
enddo
do 24 k=1,2
do 24 i=1,n(k)
x(k,i)=-aLog(amax1(.0001,ran2(idum)))/hz(k)
d(k,i)=1
temp=b+a0*ran2(idum)
if(x(k,i).gt.temp) then
x(k,i)=temp
d(k,i)=0
endif
24 continue
caLL Logrank(deLta0,test,var)
if(test.gt.za) pwr=pwr+1/fLoat(nsim)
999 continue
write(6,*)
write(6,*) '*** SimuLation resuLts (B=',nsim,')'
if(itype.eq.1) then
write(6,108) pwr,aLpha
108 format(2x,'Emp aLpha=',f7.4,' for alpha=',f6.3)
eLse
write(6,109) pwr,1-beta
109 format(2x,'Emp pwr=',f7.4,' for 1-beta=',f6.3)
endif
888 continue
stop
end
cccccccccccccccccccccccccccccccccccccccccccccccccccc
function fun(a)
common n(2),x(2,1000),d(2,1000),hz(2),p(2),
c deLta0,deLta1,b,r,za,zb,a0
externaL fv0,fv1,fw
a0=a
caLL qromb(fv0,0.,a+b,v0)
caLL qromb(fv1,0.,a+b,v1)
caLL qromb(fw,0.,a+b,w)
s0=sqrt(v0)
s1=sqrt(v1)
fun=(s0*za+s1*zb)**2/w**2-a0*r
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccc
function fv0(t)
common n(2),x(2,1000),d(2,1000),hz(2),p(2),
c deLta0,deLta1,b,r,za,zb,a0
if(t.Le.b) then
g=1
eLseif(t.ge.b.and.t.Le.a0+b) then
g=-t/a0+(a0+b)/a0
eLse
g=0
endif
s1=exp(-hz(1)*t)
s2=exp(-hz(2)*t)
temp=(p(1)*s1+deLta0*p(2)*s2)**2
fv0=g*s1*s2*(p(1)*s1*hz(1)+p(2)*s2*hz(2))/temp
fv0=fv0*deLta0*p(1)*p(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccc
function fv1(t)
common n(2),x(2,1000),d(2,1000),hz(2),p(2),
c deLta0,deLta1,b,r,za,zb,a0
if(t.Le.b) then
g=1
eLseif(t.ge.b.and.t.Le.a0+b) then
g=-t/a0+(a0+b)/a0
eLse
g=0
endif
s1=exp(-hz(1)*t)
s2=exp(-hz(2)*t)
temp=(p(1)*s1+deLta1*p(2)*s2)**2
fv1=g*s1*s2*(p(1)*s1*hz(1)+p(2)*s2*hz(2))/temp
fv1=fv1*deLta1*p(1)*p(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccc
function fw(t)
common n(2),x(2,1000),d(2,1000),hz(2),p(2),
c deLta0,deLta1,b,r,za,zb,a0
if(t.Le.b) then
g=1
eLseif(t.ge.b.and.t.Le.a0+b) then
g=-t/a0+(a0+b)/a0
eLse
g=0
endif
s1=exp(-hz(1)*t)
s2=exp(-hz(2)*t)
temp=(p(1)*s1+deLta0*p(2)*s2)*(p(1)*s1+deLta1*p(2)*s2)
fw=g*s1*s2*(p(1)*s1*hz(1)+p(2)*s2*hz(2))/temp
fw=fw*(deLta0-deLta1)*p(1)*p(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine Logrank(deLta,test,var)
reaL w(2),var,test
reaL y(2,2,10000)
common n(2),x(2,1000),d(2,1000),hz(2),p(2),
c deLta0,deLta1,b,r,za,zb,a0
do 2 k0=1,2
do 2 k=1,2
do 2 i=1,n(k)
y(k0,k,i)=0
do 2 i0=1,n(k0)
if(x(k0,i0).ge.x(k,i)) y(k0,k,i)=y(k0,k,i)+1
2 continue
var=0
do 4 k=1,2
w(k)=0
do 4 i=1,n(k)
if(d(k,i).eq.0) go to 4
w(k)=w(k)+y(3-k,k,i)/(y(1,k,i)+deLta*y(2,k,i))
var=var+y(1,k,i)*y(2,k,i)/(y(1,k,i)+deLta*y(2,k,i)) **2
4 continue
var=deLta*var
test=(deLta*w(1)-w(2))/sqrt(var)
c write(6,*) 'z=',test
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
FUNCTION ran2(idum)
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
*IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
*NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum.Le.0) then
idum=max(-idum,1)
idum2=idum
do 11 j=NTAB+8,1,-1
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.Lt.0) idum=idum+IM1
if (j.Le.NTAB) iv(j)=idum
11 continue
iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.Lt.0) idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if (idum2.Lt.0) idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.Lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
function znorm(x)
externaL erf
pi=2.*asin(1.)
z0=1.
do 2 k=1,100
f=erf(z0/sqrt(2.))/2.-.5+x
df=exp(-z0**2./2.)/sqrt(2.*pi)
znorm=z0-f/df
if(abs(znorm-z0).Lt..001.and.abs(f).Lt..001) go to 4
z0=znorm
2 continue
write(6,*) 'diverge'
4 continue
return
end
function erf(x)
if(x.Lt.0.) then
erf=-gammp(.5,x**2)
eLse
erf=gammp(.5,x**2)
endif
return
end
FUNCTION gammLn(xx)
REAL gammLn,xx
INTEGER j
REAL*8 ser,stp,tmp,x,y,cof(6)
SAVE cof,stp
DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
*24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
*-.5395239384953d-5,2.5066282746310005d0/
x=xx
y=x
tmp=x+5.5d0
tmp=(x+0.5d0)*Log(tmp)-tmp
ser=1.000000000190015d0
do 11 j=1,6
y=y+1.d0
ser=ser+cof(j)/y
11 continue
gammLn=tmp+Log(stp*ser/x)
return
END
FUNCTION gammp(a,x)
REAL a,gammp,x
CU USES gcf,gser
REAL gammcf,gamser,gLn
if(x.Lt.0..or.a.Le.0.)pause 'bad arguments in gammp'
if(x.Lt.a+1.)then
caLL gser(gamser,a,x,gLn)
gammp=gamser
eLse
caLL gcf(gammcf,a,x,gLn)
gammp=1.-gammcf
endif
return
END
SUBROUTINE gser(gamser,a,x,gLn)
INTEGER ITMAX
REAL a,gamser,gLn,x,EPS
PARAMETER (ITMAX=100,EPS=3.e-7)
CU USES gammLn
INTEGER n
REAL ap,deL,sum,gammLn
gLn=gammLn(a)
if(x.Le.0.)then
if(x.Lt.0.)pause 'x < 0 in gser'
gamser=0.
return
endif
ap=a
sum=1./a
deL=sum
do 11 n=1,ITMAX
ap=ap+1.
deL=deL*x/ap
sum=sum+deL
if(abs(deL).Lt.abs(sum)*EPS)goto 1
11 continue
pause 'a too Large, ITMAX too smaLL in gser'
1 gamser=sum*exp(-x+a*Log(x)-gLn)
return
END
SUBROUTINE gcf(gammcf,a,x,gLn)
INTEGER ITMAX
REAL a,gammcf,gLn,x,EPS,FPMIN
PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
CU USES gammLn
INTEGER i
REAL an,b,c,d,deL,h,gammLn
gLn=gammLn(a)
b=x+1.-a
c=1./FPMIN
d=1./b
h=d
do 11 i=1,ITMAX
an=-i*(i-a)
b=b+2.
d=an*d+b
if(abs(d).Lt.FPMIN)d=FPMIN
c=b+an/c
if(abs(c).Lt.FPMIN)c=FPMIN
d=1./d
deL=d*c
h=h*deL
if(abs(deL-1.).Lt.EPS)goto 1
11 continue
pause 'a too Large, ITMAX too smaLL in gcf'
1 gammcf=exp(-x+a*Log(x)-gLn)*h
return
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
SUBROUTINE qromb(func,a,b,ss)
INTEGER JMAX,JMAXP,K,KM
REAL a,b,func,ss,EPS
EXTERNAL func
PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1)
CU USES poLint,trapzd
INTEGER j
REAL dss,h(JMAXP),s(JMAXP)
h(1)=1.
do 11 j=1,JMAX
caLL trapzd(func,a,b,s(j),j)
if (j.ge.K) then
caLL poLint(h(j-KM),s(j-KM),K,0.,ss,dss)
if (abs(dss).Le.EPS*abs(ss)) return
endif
s(j+1)=s(j)
h(j+1)=0.25*h(j)
11 continue
pause 'too many steps in qromb'
END
SUBROUTINE trapzd(func,a,b,s,n)
INTEGER n
REAL a,b,s,func
EXTERNAL func
INTEGER it,j
REAL deL,sum,tnm,x
if (n.eq.1) then
s=0.5*(b-a)*(func(a)+func(b))
eLse
it=2**(n-2)
tnm=it
deL=(b-a)/tnm
x=a+0.5*deL
sum=0.
do 11 j=1,it
sum=sum+func(x)
x=x+deL
11 continue
s=0.5*(s+(b-a)*sum/tnm)
endif
return
END
SUBROUTINE poLint(xa,ya,n,x,y,dy)
INTEGER n,NMAX
REAL dy,x,y,xa(n),ya(n)
PARAMETER (NMAX=10)
INTEGER i,m,ns
REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
ns=1
dif=abs(x-xa(1))
do 11 i=1,n
dift=abs(x-xa(i))
if (dift.Lt.dif) then
ns=i
dif=dift
endif
c(i)=ya(i)
d(i)=ya(i)
11 continue
y=ya(ns)
ns=ns-1
do 13 m=1,n-1
do 12 i=1,n-m
ho=xa(i)-x
hp=xa(i+m)-x
w=c(i+1)-d(i)
den=ho-hp
if(den.eq.0.)pause 'faiLure in poLint'
den=w/den
d(i)=hp*den
c(i)=ho*den
12 continue
if (2*ns.Lt.n-m)then
dy=c(ns+1)
eLse
dy=d(ns)
ns=ns-1
endif
y=y+dy
13 continue
return
END
Fortran Ends ===================>