/*******************
This SAS program generates single-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, and beta2.
The current setting generates the example as described by the end of Section 2 in the paper. 
Users should modify the input parameters according to the actual setting.

The output parameters are:
flag =1 or 0 (1 means a success and 0 means a failure for running the program)
b (decision boundary. We used the letter "a" rather than the letter "b" in the paper)
d (decision boundary)
n (sample size required 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 Yi be the # of responders among the n patients in arm i, i=1,2. If Yi>=b then 
Arm i passes the efficacy test. If Y1-Y2>=d then Arm 1 is the winner and Arm 2 is the loser.
*******************/

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

proc iml;

start setup;
  p0_star=0.05;
  p1_star=0.2;
 * epsilon=0.000001;

  delta1=0.1;
  delta2=0.2;

  alpha1=0.05;
  alpha2=0.1;
  beta1=0.2;
  beta2=0.2;

  uu=j(10,1,0);
finish setup;

start cd(x,p,n);
/******************
P(Y<=x|p)
******************/
  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);
/*****************
P(Y=x|p)
******************/
  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 grp1_pass(b,d,n,p1,p2);
/****************************
Probability that Arm 1 passes and Arm 2 fails.
P(Y1>=b,Y2<b)+P(Y1>=b,Y2>=b,Y1-Y2>=d) 
****************************/
  y=(1-cd(b-1,p1,n))*cd(b-1,p2,n);
  do i=b to n;
    if i-d>=b then 
	  y=y+pd(i,p1,n)*(cd(i-d,p2,n)-cd(b-1,p2,n));
  end;
  return (y);
finish grp1_pass;

start grp2_pass(b,d,n,p1,p2);
/****************************
Probability that Arm 2 passes and Arm 1 fails.
P(Y2>=b,Y1<b)+P(Y2>=b,Y1>=b,Y2-Y1>=d)
****************************/
  y=grp1_pass(b,d,n,p2,p1);
  return (y);
finish grp2_pass;

start type1b(b,d,n,delta1);
/***********************
p2-p1=delta1.
Maximum of the probability that Arm 1 passes and Arm 2 fails.
*************************/
  tt=j(2,1,0);
  z=0;
  do i=1 to 99;
    p1=i*0.01;
	p2=p1+delta1;
	if p2>1 then p2=1;
	y=grp1_pass(b,d,n,p1,p2);
	if y>z then do;
      z=y;
      tt[1]=p2;
	end;
  end;
  tt[2]=z;
  return (tt);
finish type1b;

start powerb(b,d,n,p1_star,delta2);
/***********************
p2-p1=delta2 and p2>=p1_star.
Minimum of the probability that Arm 2 passes and Arm 1 fails.
*************************/
  tt=j(2,1,0);
  z=1;
  l=(1-p1_star)/100;
  do i=0 to 99;
    p2=p1_star+i*l;
	p1=p2-delta2;
	if p1<0 then p1=0;
    y=grp2_pass(b,d,n,p1,p2);
	if y<z then do;
      z=y;
      tt[1]=p2;
	end;
  end;
  tt[2]=z;
  return (tt);
finish powerb;

start search;
/**************************
Search the design that satisfies the alpha1, alpha2, beta1, and beta2 requirements.
**************************/
uu=j(10,1,0);
  al2=j(2,1,0);
  pw2=j(2,1,0);
  flag=0;
  do n=10 to 100 while (flag=0);
  do b=1 to n while (flag=0);
  do d=1 to n while (flag=0);
    al1=1-cd(b-1,p0_star,n);
	if al1<=alpha1 then do;
	  pw1=1-cd (b-1, p1_star,n);
	  if pw1>1-beta1 then do;
	    al2=type1b(b,d,n,delta1);
		if al2[2]<alpha2 then do;
		  pw2=powerb(b,d,n,p1_star,delta2);
		  if pw2[2]>=1-beta2 then flag=1; 
		  if flag=1 then do;
		    uu[1]=flag;
            uu[2]=b;
            uu[3]=d;
            uu[4]=n;
            uu[5]=al1;
            uu[6]=al2[1];
            uu[7]=al2[2];
            uu[8]=pw1;
            uu[9]=pw2[1];
			uu[10]=pw2[2];
		  end;
		end;
	  end;
	end;
  end;
  end;

  end;
finish search;

start printopt(uu);
  flag=uu[1];
     b=uu[2];
     d=uu[3];
     n=uu[4];
   al1=uu[5];
*p2_al2=uu[6];
   al2=uu[7];
   pw1=uu[8];
*p2_pw2=uu[9];
   pw2=uu[10];

print flag b d n al1 pw1 al2 pw2;

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;