/*******************
This SAS program generates two-stage designs as in the paper "Designs for Randomized 
Phase II Clinical Trials with Two Treatment Arms" by Wei Hou, myron Chang, Sin-Ho Jung, 
and Yang Li.

The input parameters are set up in the "setup" section:
           p0_star, p1_star, delta1, delta2, alpha1, alpha2, beta1, beta2, and nn,
where nn is half of the sample size needed for the corresponding single-stage design.
 
The current setting generates the example in Section 3.4 in the paper. 
Users should modify the input parameters according to the actual setting. to obtain nn, 
the user should run the single-stage design first. 

The output parameters are:
a, b, d1, and d2 (decision boundary. We used letters "a1" and "a" rather than the letters "a" and "b" in the paper)
n1 and n2 (sample sizes required for stages 1 and 2 under the significance level and power constrains)
al1 (actual type I1 error rate)
pw1 (actual power1)
al2 (actual type I2 error rate)
pw2 (actual power2)

Let Yij be the # of responders in Arm i and at stage j, i=1,2 and j=1,2. 
Y11>=a: Arm 1 passes the futility test at stage 1.
Y11-Y21>=d1: Arm 1 is the winner at stage 1.
Y11+Y12>=b: Arm 1 passes the efficacy test at stage 2.
(Y11+Y12)-(Y21+Y22)>=d2: Arm 1 is the winner at stage 2.
*******************/

options ls=80 ps=55 source center;
ods html newfile=proc;

proc iml;

start setup;
  p0_star=0.05;
  p1_star=0.2;
  nn=29;

  delta1=0.1;
  delta2=0.2;

  alpha1=0.05;
  alpha2=0.1;
  beta1=0.2;
  beta2=0.2;
  uu=j(12,1,0);

  nn=14;

finish setup;

start cd(x,p,n);
  if x<0 then y=0;
  else if x>n then y=1;
  else  do;
    u=floor(x);
	y=cdf('binomial',u,p,n);
  end;
  return (y);
finish cd;

start pd(x,p,n);
  if x<0 then y=0;
  else if x>n then y=0;
  else if x ^= floor(x) then y=0;
  else y=pdf('binomial',x,p,n);
  return (y);
finish pd;

start p_rej(a,b,n1,n2,p);
/**************************
The probability that rejecting H0 in an arm, i.e. 
P(Y1>=a, Y1+y2>=b) or 1-(P(Y1<a)+P(Y1>=a,Y1+Y2<b)). 
***************************/
  y=cd(a-1,p,n1);
  do i=a to n1;
    if b-i>=0 then y=y+pd(i,p,n1)*cd(b-i-1, p, n2);
  end;
  z=1-y;
  return (z);
finish p_rej;

start a_range(n1,p2) global(beta1);
  flag=0;
  z=0;
  do i=0 to n1+1 until (flag=1);
    y=cd(i-1,p2,n1);
	if y> beta1 then do;
	  flag=1;
	  z=i-2;
	end;
  end;
  return (z);
finish a_range;

start b_range(a,n1,n2,p1,p2) global(alpha1,beta1);
  flag1=0;
  bb=j(4,1,0);
  do i= a to n1+n2 until (flag1=1);
    y=p_rej(a,i,n1,n2,p1);
	if y<alpha1 then do;
	  bb[1]=1;
	  bb[2]=i;
	  flag1=1;
	end;
  end;

  if bb[1]=1 then do;
    u=bb[2];
	flag2=0;
	do i = n1+n2 to u by -1 until (flag2=1);
	  y=p_rej(a,i,n1,n2,p2);
	  if y>=1-beta1 then do;
	    bb[3]=1;
		bb[4]=i;
		flag2=1;
	  end;
	end;
  end;

  aa=j(3,1,0);
  if bb[1]=1 & bb[3]=1 then do;
    aa[1]=1;
	aa[2]=bb[2];
	aa[3]=bb[4];
  end;
  
  return(aa);
finish b_range;

start p00n_nnn(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms fail the futility test at stage 1.
*************************/
  y=cd(a-1,p1,n1)*cd(a-1,p2,n1);
  return (y);
finish p00n_nnn;

start p01n_n0n(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that Arm 1 fails the futility test and Arm 2 passes 
the futility test at stage 1, and Arm 2 fails the efficacy 
test at stage 2.
*************************/
  z=0;
  do i= a to n1;
    z=z+pd(i,p2,n1)*cd(b-i-1,p2,n2);
  end;
  y=cd(a-1,p1,n1)*z;
  return (y);
finish p01n_n0n;

start p10n_0nn(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that Arm 2 fails the futility test and Arm 1 passes 
the futility test at stage 1, and Arm 1 fails the efficacy 
test at stage 2.
*************************/
  y=p01n_n0n(a,b,d1,d2,n1,n2,p2,p1);
  return (y);
finish p10n_0nn;

start p113_00n(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and the equivalence 
test at stage 1, and both Arms fail the efficacy test at stage 2.
*************************/
  z=0;
  do i= a to n1;
    do j=a to n1;
	  if abs(j-i)<d1 then
      z=z+pd(i,p1,n1)*pd(j,p2,n1)*cd(b-i-1,p1,n2)*cd(b-j-1,p2,n2);
    end;
  end;
  return (z);
finish p113_00n;

start p111_0nn(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and Arm 1 is the winner
at stage 1, and Arm 1 fails the efficacy test at stage 2.
*************************/
  z=0;
  do i= a to n1;
    if i-d1 >= a then
      z=z+pd(i,p1,n1)*(cd(i-d1,p2,n1)-cd(a-1,p2,n1))*cd(b-i-1,p1,n2);
  end;
  return (z);
finish p111_0nn;

start p112_n0n(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and Arm 2 is the winner
at stage 1, and Arm 2 fails the efficacy test at stage 2.
*************************/
  y=p111_0nn(a,b,d1,d2,n1,n2,p2,p1);
  return (y);
finish p112_n0n;

start p10n_1nn(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that Arm 1 passes and Arm 2  fails the futility test at stage 1, 
and Arm 1 passes the efficacy test at stage 2.
*************************/
  z=0;
  do i= a to n1;
    z=z+pd(i,p1,n1)*(1-cd(b-i-1,p1,n2));
  end;
  y=cd(a-1,p2,n1)*z;
  return (y);
finish p10n_1nn;

start p111_1nn(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and Arm 1 is the winner 
at stage 1, and Arm 1 passes the efficacy test at stage 2.
*************************/
  z=0;
  do i= a to n1;
    if i-d1 >= a then do;
      z=z+pd(i,p1,n1)*(cd(i-d1,p2,n1)-cd(a-1,p2,n1))*(1-cd(b-i-1,p1,n2));
    end;
  end;
  return (z);
finish p111_1nn;

start p113_10n(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and the equivalence test 
at stage 1, and Arm 1 passes the efficacy test and Arm 2 fails the efficacy 
test at stage 2.
*************************/
  z=0;
  do i= a to n1;
    do j=a to n1;
	  if abs(j-i)<d1 then
      z=z+pd(i,p1,n1)*pd(j,p2,n1)*(1-cd(b-i-1,p1,n2))*cd(b-j-1,p2,n2);
    end;
  end;
  return (z);
finish p113_10n;

start p113_111(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and the equivalence test 
at stage 1, and both Arms pass the efficacy test and Arm 1 is the winner
at stage 2.
*************************/
  z=0;
  do i= a to n1;
    do j=a to n1;
	  if abs(j-i)<d1 then do;
	    do k=b-i to n2;
		  if i+k-j-d2 >= b-j then
      z=z+pd(i,p1,n1)*pd(j,p2,n1)*pd(k,p1,n2)*(cd(i+k-j-d2,p2,n2)-cd(b-j-1,p2,n2));
	    end;
	  end;
	end;
  end;
  return (z);
finish p113_111;

start p01n_n1n(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that Arm 2 passes and Arm 1  fails the futility test at stage 1, 
and Arm 2 passes the efficacy test at stage 2.
*************************/
  z=p10n_1nn(a,b,d1,d2,n1,n2,p2,p1);
  return (z);
finish p01n_n1n;

start p112_n1n(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and Arm 2 is the winner 
at stage 1, and Arm 2 passes the efficacy test at stage 2.
*************************/
  z=p111_1nn(a,b,d1,d2,n1,n2,p2,p1);
  return (z);
finish p112_n1n;

start p113_01n(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and the equivalence test 
at stage 1, and Arm 2 passes the efficacy test and Arm 1 fails the efficacy 
test at stage 2.
*************************/
  z=p113_10n(a,b,d1,d2,n1,n2,p2,p1);
  return (z);
finish p113_01n;

start p113_112(a,b,d1,d2,n1,n2,p1,p2);
/************************
Probability that both Arms pass the futility test and the equivalence test 
at stage 1, and both Arms pass the efficacy test and Arm 2 is the winner
at stage 2.
*************************/
  z=p113_111(a,b,d1,d2,n1,n2,p2,p1);
  return (z);
finish p113_112;

start p113_113(a,b,d1,d2,n1,n2,p1,p2);
/***************************
Probability that both Arms pass the two-stage testing procedure.
****************************/

  z=0;
  do i= a to n1;
    do j=a to n1;
	  if abs(j-i)<d1 then do;
	    do k=b-i to n2;
		  m=max(-d2+i+k-j+1,b-j);
		  if d2+i+k-j-1 >= m then
      z=z+pd(i,p1,n1)*pd(j,p2,n1)*pd(k,p1,n2)*(cd(d2+i+k-j-1,p2,n2)-cd(m-1,p2,n2));
	    end;
	  end;
	end;
  end;
  return (z);
finish p113_113;

start both_fail(a,b,d1,d2,n1,n2,p1,p2);
/***************************
Probability that both Arms fail the two-stage testing procedure.
****************************/
  y1=p00n_nnn(a,b,d1,d2,n1,n2,p1,p2);
  y2=p01n_n0n(a,b,d1,d2,n1,n2,p1,p2);
  y3=p10n_0nn(a,b,d1,d2,n1,n2,p1,p2);
  y4=p113_00n(a,b,d1,d2,n1,n2,p1,p2);
  y5=p111_0nn(a,b,d1,d2,n1,n2,p1,p2);
  y6=p112_n0n(a,b,d1,d2,n1,n2,p1,p2);
  z=y1+y2+y3+y4+y5+y6;
  return(z);
finish both_fail;

start grp1_pass(a,b,d1,d2,n1,n2,p1,p2);
/***************************
Probability that both Arm 1 pass and Arm 2 fail the two-stage 
testing procedure.
****************************/
  y7=p10n_1nn(a,b,d1,d2,n1,n2,p1,p2);
  y8=p111_1nn(a,b,d1,d2,n1,n2,p1,p2);
  y9=p113_10n(a,b,d1,d2,n1,n2,p1,p2);
  y10=p113_111(a,b,d1,d2,n1,n2,p1,p2);
  z=y7+y8+y9+y10;  
  return(z);
finish grp1_pass;

start grp2_pass(a,b,d1,d2,n1,n2,p1,p2);
/***************************
Probability that both Arm 2 pass and Arm 1 fail the two-stage 
testing procedure.
****************************/
  y11=p01n_n1n(a,b,d1,d2,n1,n2,p1,p2);
  y12=p112_n1n(a,b,d1,d2,n1,n2,p1,p2);
  y13=p113_01n(a,b,d1,d2,n1,n2,p1,p2);
  y14=p113_112(a,b,d1,d2,n1,n2,p1,p2);
  z=y11+y12+y13+y14;  
  return(z);
finish grp2_pass;

start both_pass(a,b,d1,d2,n1,n2,p1,p2);
/***************************
Probability that both Arms pass the two-stage testing procedure.
****************************/
  z=p113_113(a,b,d1,d2,n1,n2,p1,p2);
  return(z);
finish both_pass;

start type1_b(a,b,d1,d2,n1,n2,delta1);
/***************************
Sup of Probabilities that Arm 1 passes and Arm 2 fails the two-stage 
testing procedure when p2 -p1 >= delta1.
****************************/
  z=-1;
  do i=1 to 9;
    p1=i*0.1;
	p2=p1+delta1;
	if p2>1 then p2=1;
	y=grp1_pass(a,b,d1,d2,n1,n2,p1,p2);
	if y>z then do;
      z=y;
	  pp1=p1;
	end;
  end;

  do i=0 to 10;
    p1=pp1+(i-5)*0.01;
	p2=p1+delta1;
	if p2>1 then p2=1;
	if i^=5 then do;
	  y=grp1_pass(a,b,d1,d2,n1,n2,p1,p2);
	  if y>z then z=y;
	end;
  end;

  return(z);
finish type1_b;

start power_b(a,b,d1,d2,n1,n2,p1_star,delta2);
/***************************
Inf of Probabilities that Arm 2 passes and Arm 1 fails the two-stage 
testing procedure when p2 >= p1_star and p2 -p1 >= delta2.
****************************/
  z=1.1;
  l=(1-p1_star)/10;
  ll=l/10;
  
  do i=0 to 9;
    p2=p1_star+i*l;
	p1=p2-delta2;
	if p1<0 then p1=0; 
    y=grp2_pass(a,b,d1,d2,n1,n2,p1,p2);
	if y<z then do;
      z=y;
	  pp2=p2;
	end;
  end;

  do i=0 to 10;
    p2=pp2+(i-5)*ll;
	p1=p2-delta2;
	if p1<0 then p1=0;
	if i^=5 & p2>=p1_star then do;
      y=grp2_pass(a,b,d1,d2,n1,n2,p1,p2);
	  if y<z then z=y;
	end;
  end;
  return (z);
finish power_b;

start avess(a,b,d1,d2,n1,n2,p1,p2);
  y1=p00n_nnn(a,b,d1,d2,n1,n2,p1,p2);
  y3=p113_00n(a,b,d1,d2,n1,n2,p1,p2)
    +p113_10n(a,b,d1,d2,n1,n2,p1,p2)
    +p113_01n(a,b,d1,d2,n1,n2,p1,p2)
    +p113_111(a,b,d1,d2,n1,n2,p1,p2)
    +p113_112(a,b,d1,d2,n1,n2,p1,p2)
    +p113_113(a,b,d1,d2,n1,n2,p1,p2);
  y2=1-y1-y3;
  zz=y1*2*n1+y2*(2*n1+n2)+y3*2*(n1+n2);
  return (zz);
finish avess;

start search;
/**************************
Search the two-stage minimx design that satisfies the alpha1, alpha2, beta1, 
and beta2 requirements.
**************************/
  u=j(3,1,0);
  pp1=p0_star;
  pp2=p0_star;
  flag=0;
  do n1=nn-2 to nn+7 until (flag=1);
    n2=n1;
    zz=1000;
	flag_temp=0;
	do d1=1 to 7;
	  do d2=1 to 10;
	    xx=n1*p0_star;
		do a=1 to n1;
		  if abs(a-xx) <= 3 then do;
		  yy=(p0_star+p1_star)*(n1+n2)/2;
          do b=a to (n1+n2);
	        if abs(b-yy)<=3 then do;

            v1=p_rej(a,b,n1,n2,p0_star);
		    if v1<=alpha1 then do;
		    v2=p_rej(a,b,n1,n2,p1_star);
		    if v2>=1-beta1 then do;

              y=type1_b(a,b,d1,d2,n1,n2,delta1);
              if y<alpha2 then do; 
                z=power_b(a,b,d1,d2,n1,n2,p1_star,delta2);
                if z>=1-beta2 then do;
                  yy=avess(a,b,d1,d2,n1,n2,pp1,pp2);
				  if yy<zz then do;
				    flag_temp=1;
				    zz=yy;
                    l=p_rej(a,b,n1,n2,p0_star);
                    m=p_rej(a,b,n1,n2,p1_star);
                    uu[1]=1;
                    uu[2]=a; 
                    uu[3]=b; 
                    uu[4]=d1; 
                    uu[5]=d2; 
                    uu[6]=n1; 
                    uu[7]=n2; 
                    uu[8]=l; 
                    uu[9]=m; 
                    uu[10]=y; 
                    uu[11]=z;
					uu[12]=zz;
                  end; 
			    end;
			  end;
			  end;
			  end;
			  end;
			end;
		  end;
		end;
	  end;
	end;
	if flag_temp=1 then flag=1;
  end;
finish search;

start printopt(uu) global(p1_star, p0_star);
if uu[1]=1 then do;
   a=uu[2];
   b=uu[3];
   d1=uu[4];
   d2=uu[5];
   n1=uu[6];
   n2=uu[7];
   al1=uu[8];
   pw1=uu[9];
   al2=uu[10];
   pw2=uu[11];
   ave_ss=uu[12];
   print a b d1 d2 n1 n2;
   print al1 pw1 al2 pw2;
   p1=p1_star;
   p2=p1_star;
   pw_onepass=1-both_fail(a,b,d1,d2,n1,n2,p1,p2);
   p1=p0_star;
   p2=p0_star;
   al_onepass=1-both_fail(a,b,d1,d2,n1,n2,p1,p2);
   pw3_pw1=1-(1-pw1)**2-pw_onepass;
   al3_al1=2*al1-al1**2-al_onepass;
   *print pw_onepass al_onepass pw3_pw1 al3_al1 ave_ss;
 end;
 else 
   print 'no solusion';
 finish printopt;

 
start main;
  run setup;
  print alpha1 beta1 alpha2 beta2;
  print p0_star p1_star delta1 delta2;
  run search;
  run printopt(uu);
finish main;  
 
run;
quit;