ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Author: Tang Man-Lai
c Date: 7-28-1998
c Program: uncon1.for
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       double precision fact(0:50000)
       double precision cnr
       double precision alpha
       double precision ca, temp1, temp2, dbar
       double precision c(0:1,0:1000,0:10000)
       double precision rnd(10), rny(10), rnn(10)
       double precision rns, rnt, rntotal
       double precision pvmax, obsca, sdsq, ds, dt
       double precision p(0:10000), pv(0:10000)
       double precision camax, camin
       double precision pt, pv1, pb
       double precision coeff(0:500000)
       double precision gamma, beta, power
       integer ncoeff(2,0:500000)
       integer nadd(0:500000)
       integer nd(10), ny(10), nn(10)
       integer imaxs(0:1,0:200)
       integer time,n1,n2,n3
       open(unit = 8, file = 'rca.dat')
       open(unit = 16, file = 'ca.dat')

c Maximum value for the statistic S
       ismax = 1000

c Maximum value for the statistic T
       itmax = 10000

c Maximum value for factorial calculation
       imax = 50000

c The number of grid points for uniform search in identifying the
c exact unconditional p-value or critical value
       read(8,*) inc 

c iwant = 0 -> exact unconditional p-value calculation
c iwant = 1 -> exact unconditional power calculation
       read(8,*) iwant

c The nominal significance level
       read(8,*) alpha

c Number of dose levels in the study
       read(8,*) ndose

c Dose scores
       read(8,*) (nd(k), k = 1, ndose)

c Total number of observations in each dose level
       read(8,*) (nn(k), k = 1, ndose)

c If iwant= 0, then the number of successes in each dose level is
c required
       if (iwant .eq. 0) then
          read(8,*) (ny(k), k = 1, ndose)
       endif

c If iwant = 1, then the values for gamma and beta is required
       if (iwant .eq. 1) then 
          read(8,*) gamma, beta
       endif

c itest = 1 -> logistic function
c itest = 2 -> probit function 
c itest = 3 -> extreme value function
c itest = 4 -> one hit function
       if (iwant .eq. 1) then 
          read(8,*) itest
       endif
       
c
c Beginning of the main program
c

       n1 = time()

c
c Calculate preliminary statistics
c

       ntotal = 0
       ns = 0
       nt = 0
       dbar = 0.0
       do 10 k = 1, ndose
          rnd(k) = dble(dfloat(nd(k)))
          rnn(k) = dble(dfloat(nn(k)))
          ntotal = ntotal + nn(k)
          dbar = dbar + nd(k)*nn(k)
          if (iwant .eq. 0) then
             rny(k) = dble(dfloat(ny(k)))
             ns = ns + ny(k)
             nt = nt + ny(k)*nd(k)
          endif
10     continue

       nhalf = int(ntotal/2)
       rntotal = dble(dfloat(ntotal))
       dbar = dbar/rntotal
       if (iwant .eq. 0) then
          rns = dble(dfloat(ns))
          rnt = dble(dfloat(nt))
       endif

c
c Construct the c(s, t) matrix
c

       do 100 i = 0, imax
          fact(i) = 0.0
100    continue
       call factorial(fact, imax)

       io = 0
       do 110 i1 = 0, 1
          do 110 i2 = 0, nhalf
             do 110 i3 = 0, itmax
                c(i1, i2, i3) = 0.0
110    continue
       c(io, 0, 0) = 1.0
       do 111 i2 = 1, nhalf
          do 111 i1 = 0, 1
             imaxs(i1, i2) = -1
111    continue
       imaxs(io, 0) = 0

       do 120 k = 1, ndose
          in = 1 - io          
          do 130 iy = 0, nn(k)
             temp1 = cnr(nn(k), iy, fact, imax)
             ntemp2 = nd(k)*iy
             do 140 i1 = 0, nhalf-iy
                ntemp1 = i1 + iy
                do 150 i2 = 0, imaxs(io, i1)
                   if (c(io, i1, i2) .eq. 0.0) goto 150
                   c(in, ntemp1, i2+ntemp2) = 
     *             c(in, ntemp1, i2+ntemp2) + c(io, i1, i2)*temp1
                   if (i2 + ntemp2 .gt. imaxs(in, ntemp1)) then
                      imaxs(in, ntemp1) = i2 + ntemp2
                   endif 
150             continue
140          continue
130       continue
          io = in
          in = 1 - io
          do 160 i2 = 0, nhalf
             imaxs(in, i2) = -1
             do 160 i3 = 0, itmax
                c(in, i2, i3) = 0.0
160       continue
120    continue
       n2 = time()
       write(*, *) n2 - n1
       
       if (iwant .eq. 1) goto 1111

c
c (1) Determine the exact unconditional p-value
c
c
c Calculate the observed Cochran-Armitage trend test statistic
c

       temp1 = 0.0
       temp2 = 0.0
       do 20 k = 1, ndose
          temp1 = temp1 + rnd(k)*(rny(k) - rns*rnn(k)/rntotal)
          temp2 = temp2 + rnn(k)*((rnd(k) - dbar)**2.0)
20     continue
       temp1 = temp1**2.0
       sdsq = temp2
       temp2 = temp2*rns*(rntotal - rns)/rntotal/(rntotal - 1.0)
       if (temp2 .ne. 0.0) then
          ca = temp1/temp2
       else
          ca = 0.0
       endif
c       write(*,*) ca

       do 300 j = 0, inc       
          p(j) = dble(dfloat(j))/dble(dfloat(inc))/2.0
          pv(j) = 0.0
300    continue

       do 310 i1 = 0, nhalf
          ds = dble(dfloat(i1))
          do 320 i2 = 0, imaxs(io, i1)   
             if (c(io, i1, i2) .eq. 0.0) goto 320
             dt = dble(dfloat(i2))
             obsca = (dt - ds*dbar)**2.0
             if (obsca .ne. 0.0) then
                obsca = (obsca*rntotal*(rntotal - 1.0))/sdsq/ds/
     *                  (rntotal - ds)
             endif
             if (obsca .ge. ca) then
                do 330 j = 0, inc
                   if ((i1 .eq. nhalf) .and. (2*nhalf .eq. ntotal)) then
                      pv(j) = pv(j) + c(io, i1, i2)*
     *                     (p(j)**ds)*((1.0 - p(j))**(rntotal - ds)) 
                   else
                      pv(j) = pv(j) + c(io, i1, i2)*
     *                     ((p(j)**ds)*((1.0 - p(j))**(rntotal - ds)) +
     *                      (p(j)**(rntotal - ds))*((1.0 - p(j))**ds))
                   endif
330             continue
             endif
320       continue
310    continue

       pvmax = 0.0
       do 340 j = 0, inc
          if (pv(j) .gt. pvmax) pvmax = pv(j)
340    continue
       write(*,*) pvmax

       n3 = time()
       write(*, *) n3 - n2
       goto 999

c
c (2) Determine the exact unconditional power
c

c
c Sort the Cochran-Armitage trend statistic values
c

1111   temp2 = 0.0
       do 400 k = 1, ndose
          temp2 = temp2 + rnn(k)*((rnd(k) - dbar)**2.0)
400     continue
       sdsq = temp2

       icount = 0
       do 410 i1 = 0, nhalf
          ds = dble(dfloat(i1))
          do 420 i2 = 0, imaxs(io, i1)   
             if (c(io, i1, i2) .eq. 0.0) goto 420
             dt = dble(dfloat(i2))
             temp1 = (dt - ds*dbar)**2.0
             if (temp1 .ne. 0.0) then
                temp1 = (temp1*rntotal*(rntotal - 1.0))/sdsq/ds/
     *                  (rntotal - ds)
             endif
             if (icount .eq. 0) then
                icount = icount + 1
                coeff(icount) = temp1
                ncoeff(1,icount) = i1
                ncoeff(2,icount) = i2
                nadd(icount) = -999
                nmax = icount
                camax = temp1
                nmin = icount
                camin = temp1
             else
                if (temp1 .ge. camax) then
                   icount = icount + 1
                   coeff(icount) = temp1
                   ncoeff(1,icount) = i1
                   ncoeff(2,icount) = i2
                   nadd(icount) = nmax
                   nmax = icount
                   camax = temp1
                else if (temp1 .le. camin) then
                   icount = icount + 1
                   coeff(icount) = temp1
                   ncoeff(1,icount) = i1
                   ncoeff(2,icount) = i2
                   nadd(nmin) = icount 
                   nmin = icount
                   camin = temp1
                   nadd(icount) = -999
                else
                   nbefore = nmax
                   nafter = nadd(nbefore)
444                if ((temp1 .le. coeff(nbefore)) .and.
     *                 (temp1 .ge. coeff(nafter))) then
                      icount = icount + 1
                      coeff(icount) = temp1
                      ncoeff(1,icount) = i1
                      ncoeff(2,icount) = i2
                      nadd(nbefore) = icount
                      nadd(icount) = nafter
                   else
                      nbefore = nafter
                      nafter = nadd(nafter)
                      goto 444
                   endif
                endif    
             endif
420       continue
410    continue

c
c Identify the critical value 
c

       do 500 j = 1, inc       
          pt = dble(dfloat(j))/dble(dfloat(inc))/2.0
          pv1 = 0.0
          ncurrent = nmax
          pb = coeff(ncurrent)
 555      i1 = ncoeff(1,ncurrent)
          i2 = ncoeff(2,ncurrent)
          ds = dble(dfloat(i1))
          dt = dble(dfloat(i2))
          if ((i1 .eq. nhalf) .and. (2*nhalf .eq. ntotal)) then
             pv1 = pv1 + c(io, i1, i2)*
     *                     (pt**ds)*((1.0 - pt)**(rntotal - ds)) 
          else
             pv1 = pv1 + c(io, i1, i2)*
     *                     ((pt**ds)*((1.0 - pt)**(rntotal - ds)) +
     *                      (pt**(rntotal - ds))*((1.0 - pt)**ds))
          endif
          ntemp1 = nadd(ncurrent)
          if (ntemp1 .eq. -999) goto 556
          if (coeff(ntemp1) .eq. coeff(ncurrent)) then
             ncurrent = ntemp1
             goto 555
          endif

556       if (pv1 .ge. alpha) then
             pv(j) = pb
             goto 500
          endif

          pb = coeff(ncurrent)
          ncurrent = ntemp1
          if (ncurrent .ne. -999) goto 555
500    continue
   
       pvmax = 0
       do 510 j = 0, inc
          if (pv(j) .gt. pvmax) pvmax = pv(j)
510    continue

       if (itest .gt. 1) goto 777
c
c Calculate the exact unconditional power based on linear logistic model
c
       temp1 = 1.0
       do 600 k = 1, ndose
          p(k) = dexp(gamma + beta*rnd(k))/
     *           (1.0 + dexp(gamma + beta*rnd(k)))
          temp1 = temp1*(1.0 + dexp(gamma + beta*rnd(k)))**rnn(k)
          write(16,*) p(k)
600    continue

       power = 0.0
       ncurrent = nmax
666    if (coeff(ncurrent) .ge. pvmax) then
          i1 = ncoeff(1,ncurrent)
          i2 = ncoeff(2,ncurrent)
          ds = dble(dfloat(i1))
          dt = dble(dfloat(i2))
          if ((i1 .eq. nhalf) .and. (2*nhalf .eq. ntotal)) then
             power = power + c(io, i1, i2)*
     *               dexp(gamma*ds + beta*dt) 
          else
             power = power + c(io, i1, i2)*
     *               (dexp(gamma*ds + beta*dt) +
     *               dexp(gamma*(rntotal - ds) + 
     *                    beta*(rntotal*dbar - dt)))
          endif
          ncurrent = nadd(ncurrent)
          goto 666
       endif
       power = power/temp1
       write(*,*) power

       n3 = time()
       write(*, *) n3 - n2

       goto 999

c
c Calculate the exact unconditional power based on other linear models
c

777    power = 0.0
       do 700 k = 1, ndose
          temp1 = gamma + beta*rnd(k)
          if (itest .eq. 2) then
             call dmdnor(temp1,p(k))
          endif
          if (itest .eq. 3) then
             p(k) = 1.0 - dexp(-1.0*dexp(temp1))
          endif
          if (itest .eq. 4) then
             p(k) = 1.0 - dexp(-1.0*temp1)
          endif
          if (itest .eq. 5) then
             p(k) = dexp(temp1)/(1.0 + dexp(temp1))
          endif
          write(16,*) p(k)
700    continue

c
c Construct the first half of the probability generating function
c

       io = 0
       do 710 i1 = 0, 1
          do 710 i2 = 0, nhalf
             do 710 i3 = 0, itmax
                c(i1, i2, i3) = 0.0
710    continue
       c(io, 0, 0) = 1.0
       do 711 i2 = 1, nhalf
          do 711 i1 = 0, 1
             imaxs(i1, i2) = -1
711    continue
       imaxs(io, 0) = 0

       do 720 k = 1, ndose
          in = 1 - io          
          do 730 iy = 0, nn(k)
             temp1 = cnr(nn(k), iy, fact, imax)*
     *               (p(k)**iy)*((1.0 - p(k))**(nn(k) - iy))
             ntemp2 = nd(k)*iy
             do 740 i1 = 0, nhalf-iy
                ntemp1 = i1 + iy
                do 750 i2 = 0, imaxs(io, i1)
                   if (c(io, i1, i2) .eq. 0.0) goto 750
                   c(in, ntemp1, i2+ntemp2) = 
     *             c(in, ntemp1, i2+ntemp2) + c(io, i1, i2)*temp1
                   if (i2 + ntemp2 .gt. imaxs(in, ntemp1)) then
                      imaxs(in, ntemp1) = i2 + ntemp2
                   endif 
750             continue
740          continue
730       continue
          io = in
          in = 1 - io
          do 760 i2 = 0, nhalf
             imaxs(in, i2) = -1
             do 760 i3 = 0, itmax
                c(in, i2, i3) = 0.0
760       continue
720    continue

       ncurrent = nmax
770    if (coeff(ncurrent) .ge. pvmax) then
          i1 = ncoeff(1,ncurrent)
          i2 = ncoeff(2,ncurrent)
          power = power + c(io, i1, i2)
          ncurrent = nadd(ncurrent)
          goto 770
       endif
       
c
c Construct the second half of the probability generating function
c

       io = 0
       do 810 i1 = 0, 1
          do 810 i2 = 0, nhalf
             do 810 i3 = 0, itmax
                c(i1, i2, i3) = 0.0
810    continue
       c(io, 0, 0) = 1.0
       do 811 i2 = 1, nhalf
          do 811 i1 = 0, 1
             imaxs(i1, i2) = -1
811    continue
       imaxs(io, 0) = 0

       do 820 k = 1, ndose
          in = 1 - io          
          do 830 iy = 0, nn(k)
             temp1 = cnr(nn(k), iy, fact, imax)*
     *               ((1.0 - p(k))**iy)*(p(k)**(nn(k) - iy))
             ntemp2 = nd(k)*iy
             do 840 i1 = 0, nhalf-iy
                ntemp1 = i1 + iy
                do 850 i2 = 0, imaxs(io, i1)
                   if (c(io, i1, i2) .eq. 0.0) goto 850
                   c(in, ntemp1, i2+ntemp2) = 
     *             c(in, ntemp1, i2+ntemp2) + c(io, i1, i2)*temp1
                   if (i2 + ntemp2 .gt. imaxs(in, ntemp1)) then
                      imaxs(in, ntemp1) = i2 + ntemp2
                   endif 
850             continue
840          continue
830       continue
          io = in
          in = 1 - io
          do 860 i2 = 0, nhalf
             imaxs(in, i2) = -1
             do 860 i3 = 0, itmax
                c(in, i2, i3) = 0.0
860       continue
820    continue

       ncurrent = nmax
870    if (coeff(ncurrent) .ge. pvmax) then
          i1 = ncoeff(1,ncurrent)
          i2 = ncoeff(2,ncurrent)
          if ((i1 .ne. nhalf) .or. 
     *        ((i1 .eq. nhalf) .and. (2*nhalf .ne. ntotal))) then
              power = power + c(io, i1, i2)
          endif
          ncurrent = nadd(ncurrent)
          goto 870
       endif
       write(*,*) power

       n3 = time()
       write(*, *) n3 - n2

999    stop
       end

ccccccccccccccccccccccccccccccccccc
c
c Calculate the Factorial!
c
ccccccccccccccccccccccccccccccccccc
       subroutine factorial(fact, imax)
       double precision fact(0:imax)
       fact(0) = 0.0
       do 100 i = 1, imax
          fact(i) = fact(i-1) + dlog(dble(float(i)))
100    continue
       return
       end

cccccccccccccccccccccccccccccccccccc
c
c Calculate C(n,r)!
c
cccccccccccccccccccccccccccccccccccc
       double precision function cnr(n,nr,fact,imax)
       double precision fact(0:imax)
       n_nr = n - nr
       cnr = fact(n)-fact(nr)-fact(n_nr)
       cnr = dexp(cnr)
       return
       end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c function.dmdnor
c
c this is to compute the prob. that the r.v. x less than the value
c y where x is distributed as a standard normal distribution
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine dmdnor(y,p)
      double precision y,p,sqr1d2,erfcc
      sqr1d2 = 0.70710678118654d0
      p = -y*sqr1d2
      if (dabs(p) .le. 13.2d0) goto 5
      p = 0.0d0
      if (y .lt. 0.0d0) return
      p = 1.0d0
      return
5     p = 0.5*erfcc(p)
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c This function returns the complementary error function
c erfc(x) with fractional error everywhere less than
c 1.2d-07.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision function erfcc(x)
      double precision x,z,t
      z = dabs(x)
      t = 1.0/(1.0+0.5*z)
      erfcc = t*dexp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
     *      t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
     *      t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
      if (x .lt. 0.00) erfcc=2.0-erfcc
      return
      end
 
 

100
1
0.05
3
0    1    2       
5   5   5  
-3.0    0.1
2