/*Two-stage design for a phase II window clinical trial.sas.

  This program generates designs for phase II window trails as described in 
  the paper entitled "One- and two-stage designs for phase II window studies",
  Statistics in Medicine, 2007, 26:2604-2614, by Chang, MN, M devidas 
  and James Andersons
 
  The hypotheses are H0: p1<=p01 or p2>=p02 vs. H1: p1>=p11 and p2<=p12, 
  where p01<p11 and p02>p12. 

  The rejection region of H0 for a single stage test is that (S1>=s1 and T1<=t1)
  or (continue to the second stage, and S>=a1 and T<=a2).

  The input parameters are p01, p11, p02, p12, alpha, and beta (1 - power) 
  as describe in the paper. In addiiton, aalpha is Z-sub-alpha and zbeta is 
  Z-sub_beta. For example, if alpha = 0.05 and beta = 0.2, then zalpha = 1.645 
  and zbeta=0.84. All these input parameters should be specified in the setup 
  section by the user. The current setting is corresprding to the first row 
  entry of Table II in the paper.
  
  The user should run this program without change of the input parameters first 
  and compare the results with the output parameters in the first row entry of 
  Table II in the paper. By this way the user will understand the meaning of the 
  output parameters generated by the SAS program. 
*/

options ls=80 ps=55 source center;

proc iml;

start setup;
  p01=0.1;
  p11=0.3;
  p02=0.9;
  p12=0.7;
  alpha=0.05;
  beta=0.2;
  zalpha=1.645;
  zbeta=0.84;

  pdf1=j(4,100,0);
  cdf1=j(4,100,1);
  cdf2=j(4,100,1);
  cc=j(13,1,-1);

  cnst=0.009;
  neg=10**(-10);
  one=1+2*neg;
finish setup;

start power1(m,s,t) global(p11,p12,neg);
/*
  This is the power function of rejecton region S>=s and T<=t for 
  a single-stage test with response rate of p1 and early progression 
  rate of p2. Constrains: 0<p1<1; 0<p2<1; p1+p2<=1; s+t<=m.
*/ 

  if ((t<0) | (s>m)) then y=0;
  else if (t>m) then y=1-cdf('binomial',s-1,p11,m);
  else if (s<0) then y=cdf('binomial',t,p12,m);
  else do; 
    y=0;
    do i=0 to t;
      y=y+pdf('binomial',i,p12,m)*(1-cdf('binomial',s-1,p11/(1-p12)-neg,m-i));
    end;
  end;

return (y);
finish power1;

start mul(n,s,t,p1,p2) global(neg);
* This is the multinomial probability;
*print n s t p1 p2;
  if (s<0) | (t<0) | (s+t>n) then y=0;
  else
    y=pdf('binomial',s,p1,n)*pdf('binomial',t,p2/(1-p1)-neg,n-s);

return(y);
finish mul;

start contn(n1,s1,s2,t1,t2,s,t);
/*
  contn=1 means that (s,t) belongs to the continuation region.
  The rejection region at the first stage is that S1>=s1 and 
  T1<=t1; The acceptance region is that S1<=s2 or T1>=t2.
  It is assumed that s1>s2 and t1<t2.
*/

  y=1;
  if (s>=s1) & (t<=t1) then y=0;
  else if (s<=s2) | (t>=t2) then y=0;
  else if (s+t>n1) | (s<0) | (t<0) then y=0;

return (y);
finish contn;


start power2(n1,n2,s1,s2,t1,t2,a1,a2) global(p11,p12,neg);
/*
  This subroutine computes the probability that (S1>=s1 and T1<=t1) or 
  {[(T1<t2 and s2<S1<s1) or (t1<T1<t2 and s1<=S1<=n1-T1)]
  and [S>=a1 and T<=a2]}.
*/

  y=power1(n1,s1,t1);
  do s=0 to n1;
    do t=0 to n1-s;
      flag=contn(n1,s1,s2,t1,t2,s,t);
      if (flag=1) then
        y=y+mul(n1,s,t,p11,p12)*power1(n2,a1-s,a2-t);
    end;
  end;

return(y);
finish power2;


start sz(p0,p1) global(alpha, beta);
/*
  This subroutine computes the exact sample size for a binomial test:
  H0: p=p0 vs. H1: p=p1
*/

  flag=0;
  m=1;
  do until (flag=1);
    al=0;
    s=m+1;
    do until ((al>alpha) | (s=0));
       s=s-1;
       al=1-cdf('binomial',s-1,p0,m);
    end;
    s=s+1;

    bet=cdf('binomial',s-1,p1,m);
    if bet<=beta then do;
      flag=1;
      x=m;
    end;
    else m=m+1;
  end;

return (x);
finish sz;

start power_reg(a,b,c,dd) global(pdf1,cdf1,cdf2);
/*
The testing procedure is to reject H0 if (S1>=b) or ((a+1<=S1<=b-1) 
and (S1+S2>=c)).

If dd=1, then the subroutine computes the power of the test specified 
              by the stopping boundaries a,b, and c under p=p01;
If dd=2, then the subroutine computes the power of the test specified 
              by the stopping boundaries a,b, and c under p=p11;
If dd=3, then the subroutine computes the power of the test specified 
              by the stopping boundaries a,b, and c under p=1-p02;
If dd=4, then the subroutine computes the power of the test specified 
              by the stopping boundaries a,b, and c under p=1-p12;

For definition of pdf1[ ], cdf1[ ], and cdf2[ ], see subroutine "search". 

*/
 
  pdff=j(100,1,0);
  cdff1=j(100,1,1);
  cdff2=j(100,1,1);

  if dd=1 then do;
    pdff=pdf1[1,];
    cdff1=cdf1[1,];
    cdff2=cdf2[1,];
  end;
  if dd=2 then do;
    pdff=pdf1[2,];
    cdff1=cdf1[2,];
    cdff2=cdf2[2,];
  end;
  if dd=3 then do;
    pdff=pdf1[3,];
    cdff1=cdf1[3,];
    cdff2=cdf2[3,];
  end;
  if dd=4 then do;
    pdff=pdf1[4,];
    cdff1=cdf1[4,];
    cdff2=cdf2[4,];
  end;

  if b>0 then
    y=1-cdff1[b];
  else y=1;

  do i=a+1 to b-1;
    if c-i>0 then 
      x=1-cdff2[c-i];
    else x=1;
    y=y+pdff[i+1]*x;
  end;

  return (y);
  finish power_reg;


start eval(ee) global(n1,n2,alpha,beta,pdf1,cdf1,cdf2,p01,p02,p11,p12,cnst);
/*
If ee=1, then eval(1) is the set of designs with sample sizes n1 and n2 
              satisfying the sig and power requirement for H0: p=p01 vs H1: p=p11;
If ee=2, then eval(2) is the set of designs with sample sizes n1 and n2 
              satisfying the sig and power requirement for H0: p=1-p02 vs H1: p=1-p12;

For definition of cdf1[ ], see subroutine "search".

*/
 
  cdff1=j(100,1,1);
  cdff2=j(100,1,1);

  if ee=1 then do;
    cdff0=cdf1[1,];
    cdff1=cdf1[2,];
  end;
  if ee=2 then do;
    cdff0=cdf1[3,];
    cdff1=cdf1[4,];
  end;

  x=j(6,400,0);
  dd0=1+(ee-1)*2;
  dd1=2+(ee-1)*2;

  ngl=min(alpha,beta)*cnst;

  u1=-1;
  bet1=0;
  do until ((bet1>ngl) | (u1=n1));
    u1=u1+1;
    bet1=min(cdff1[u1+1],cdff0[u1+1]);
  end;
  if bet1>ngl then u1=u1-1;
  else u1=n1;

  u2=n1+1;
  al1=0;
  do until ((al1>ngl) | (u2=1));
    u2=u2-1;
    al1=min(1-cdff0[u2],1-cdff1[u2]);
  end;
  if al1>ngl then u2=u2+1;
  else u2=1;

/*
  To save computing time, [0,u1-1] and [u2+1,n1] will not be considered 
  as candidate for a and b, respectively, since the prob is negligible.
*/

  m1=-1;
  bet1=0;
  do until ((bet1>beta) | (m1=n1));
    m1=m1+1;
    bet1=cdff1[m1+1];
  end;
  if bet1>beta then m1=m1-1;
  else m1=n1;

  m2=n1+1;
  al1=0;
  do until ((al1>alpha) | (m2=1));
    m2=m2-1;
    al1=1-cdff0[m2];
  end;
  if al1>alpha then m2=m2+1;
  else m2=1;

*print m1 u1 m2 u2;

  if m2<=m1 then do;
    x[1,1]=-1;
    goto fineval;
/* 
   If m2<=m1, then the single stage test with sample size n1 and critical 
   value=(m1+m2)/2 will satisfy the power and sig requirements.
*/
  end;
  else do;
  
  k=0;
  do i=m1 to u1 by -1;
  do j=m2 to u2;
    if i<j then do;
      a=0;
      b=n1+n2;
      do until (b-a<=1);
        mid=int((a+b)/2);
        u=power_reg(i,j,mid,dd0);
        v=power_reg(i,j,mid,dd1);
*print i j mid u v;
        if u>=alpha & v<=1-beta then goto contloop;
        else if u>=alpha & v>1-beta then a=mid;
        else if u<alpha & v<=1-beta then b=mid;
        else if u<alpha & v>1-beta & k<400 then do;

u0=u;
do until ((u0>alpha) | (mid=a));
  mid=mid-1;
  u0=power_reg(i,j,mid,dd0);
end;
mid=mid+1;

          k=k+1;
          if ee=1 then do;
            x[1,1]=k;
            x[2,k]=i;
            x[3,k]=j;
            x[4,k]=mid;
            x[5,k]=v;
            x[6,k]=u;
          end;
          else if ee=2 then do;
            x[1,1]=k;
            x[2,k]=n1-j;
            x[3,k]=n1-i;
            x[4,k]=n1+n2-mid;
            x[5,k]=v;
            x[6,k]=u;
          end;
          goto contloop;
        end;
      end;
    end;
contloop:
  end;
  end;
  end;

fineval:
  return (x);
  finish eval;

start avesz(n1,n2,s1,s2,t1,t2,p1,p2) global(neg);
/*
  This subroutine computes the average sample size under the parameters
  (p1,p2) and with the stopping boundaries s1, s2, t1, and t2 at the 
  first stage. 
*/

  y=n1;
  do s=0 to n1;
    do t=0 to n1-s;
      z=contn(n1,s1,s2,t1,t2,s,t);
      if z=1 then y=y+mul(n1,s,t,p1,p2)*n2;
    end;
  end;
  return (y);
finish avesz;


start search;
/*
  This subroutine finds the sample sizes and stopping boundaries
  for a 2_stage test with the minimum ave. sample size.

  The rejection region of H0 is: (S1>=s1 and T1<=t1) or 
  {[(T1<t2 and s2<S1<s1) or (t1<T1<t2 and s1<=S1<=n1-T1)]
  and [S>=a1 and T<=a2]}.

  pdf1[1, ]: binomial probs with n1 and p01;
  pdf1[2, ]: binomial probs with n1 and p11;
  pdf1[3, ]: binomial probs with n1 and 1-p02;
  pdf1[4, ]: binomial probs with n1 and 1-p12;
  cdf1[1, ]: binomial distribution with n1 and p01;
  cdf1[2, ]: binomial distribution with n1 and p11;
  cdf1[3, ]: binomial distribution with n1 and 1-p02;
  cdf1[4, ]: binomial distribution with n1 and 1-p12;
  cdf2[1, ]: binomial distribution with n2 and p01;
  cdf2[2, ]: binomial distribution with n2 and p11;
  cdf2[3, ]: binomial distribution with n2 and 1-p02;
  cdf2[4, ]: binomial distribution with n2 and 1-p12;

*/

  m1=sz(p01,p11);
  m2=sz(1-p02,1-p12);
*print m1 m2;

  m3=max(m1,m2);
  n1=int(m3/2);
  if n1>3 then 
    n2=n1-3;
  else n2=n1;

  flag=0;
  do until ((flag=1) | (n1=50));

*print n1 n2;

  do i=0 to n1;
    pdf1[1,i+1]=pdf('binomial',i,p01,n1);
    pdf1[2,i+1]=pdf('binomial',i,p11,n1);
    pdf1[3,i+1]=pdf('binomial',i,1-p02,n1);
    pdf1[4,i+1]=pdf('binomial',i,1-p12,n1);
    cdf1[1,i+1]=cdf('binomial',i,p01,n1);
    cdf1[2,i+1]=cdf('binomial',i,p11,n1);
    cdf1[3,i+1]=cdf('binomial',i,1-p02,n1);
    cdf1[4,i+1]=cdf('binomial',i,1-p12,n1);
  end;
  do i=0 to n2;
    cdf2[1,i+1]=cdf('binomial',i,p01,n2);
    cdf2[2,i+1]=cdf('binomial',i,p11,n2);
    cdf2[3,i+1]=cdf('binomial',i,1-p02,n2);
    cdf2[4,i+1]=cdf('binomial',i,1-p12,n2);
  end;

  ave=10000;
  x=j(6,400,0);
  y=j(6,400,0);  
  x=eval(1);
*print x;
  if x[1,1]>=1 then do;
    y=eval(2);
*print y;
    if y[1,1]>=1 then do;
      m1=x[1,1];
      m2=y[1,1];
      do i=1 to m1;
        do j=1 to m2;
          s2=x[2,i];
          s1=x[3,i];
          a1=x[4,i];
          t1=y[2,j];
          t2=y[3,j];
          a2=y[4,j];
          z=power2(n1,n2,s1,s2,t1,t2,a1,a2);

          if z>=1-beta then do;
            flag=1;
            ave1=(avesz(n1,n2,s1,s2,t1,t2,p01,p02)+
                  avesz(n1,n2,s1,s2,t1,t2,p11,p12))/2;
            if ave1<ave then do;
              ave=ave1;
              cc[1]=n1;
              cc[2]=n2;
              cc[3]=s1;
              cc[4]=s2;
              cc[5]=t1;
              cc[6]=t2;
              cc[7]=a1;
              cc[8]=a2;
              cc[9]=x[6,i];
              cc[10]=y[6,j];
              cc[11]=max(x[6,i],y[6,j]);
              cc[12]=z;
              cc[13]=ave;
            end;
          end;
        end;
      end;
    end;
  end;
    if flag=0 then do;
      if n1>n2 then n2=n2+1;
      else do;
        n1=n1+1;
        n2=n2-2;
      end;
    end;
  end;
finish search;

start main;

  run setup;
  print alpha beta;
      cc=j(13,1,-1);
      run search;
        n1=cc[1];
        n2=cc[2];
        s1=cc[3];
        s2=cc[4];
        t1=cc[5];
        t2=cc[6];
        a1=cc[7];
        a2=cc[8];
        al1=cc[9];
        al2=cc[10];
        al=cc[11];
        bet=cc[12];
        ave=cc[13];

        print p01 p11 p02 p12;
        print n1 n2 s1 s2 t1 t2 a1 a2;
        print al bet ave;

finish main;

run;