/*******************
Category_I_redesign.sas
The two-stage procedure is testing H0: p<=p0 vs. p>=p1, where p0<p1.
The numbers of patients in stages 1 and 2 are n1 and n2, respectively. n=n1+n2.
This program generates planned designs with either n1=n2 or n1=n2+1; 
The numbers of responses at stages 1 and 2 are Y1 and Y2, respectively. Y=Y1+Y2. 
If Y1<=a then the drug is rejected, and the trial is stopped at hte first stage. 
If Y1>a then the accrual is continued to the second stage. If Y<=c at the second stage 
then the drug is rejected. The significance level is: 
                           P(Y1>a, Y>c|n1,n2,p0).
The power is:              
                           P(Y1>a, Y>c|n1,n2,p1).
The average sample size is: 
                           aven=P(Y1>a|p0)*n2+n1.
The program search all possible combinaitons of (a,c,n1,n) to generate a planned 
design satisfying significance level and power requirement and minimizing the 
average sample size. 

When the actual sample sizes n1_star at stage 1 is different from the planned sample size 
n1, we generate a modified design (a_star, c_star, n1_star, n_star) for given n1_star 
satisfying the significance and power requirements and minimizing the average sample size. 
When the actual sample size n_star_star is defferent with n_sar at stage 2, we find c_star_star
such that the design(a_star, c_star_star, n1_star, n_star_star) satisfies the significance 
requirement.
*******************/

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

proc iml;

start setup;
  epsilon=0.00000001;
  p0=0.1;
  p1=0.3;
  alpha=0.05;
  beta=0.2;
finish setup;

start cd(x,p,n) global(epsilon);
/******
Compute P(Y<=x|n,p), where Y~Bin(n,p).
******/
  if p<epsilon then p=epsilon;
  if p>1-epsilon then p=1-epsilon;
  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) global(epsilon);
/******
Compute P(Y=x|n,p), where Y~Bin(n,p).
******/
  if p<epsilon then p=epsilon;
  if p>1-epsilon then p=1-epsilon;
  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 a_up(p1,n1) global(beta);
/***************
Find the largest integer uu such that
   P(Y1<=uu)<=beta
***************/
  flag=0;
  m=-1;
  do while((flag=0) & (m< n1+1));
    y=cd(m,p1,n1);
	if y<=beta then m=m+1;
	else do;
	  flag=1;
	  uu=m-1;
	end;
  end;
  return (uu);
finish a_up;
     
start power(a,c,n1,n,p);
/***********
Compute the power for given design (a,c,n1,n) under response rate p:
  P(y1>a, Y>c|p,n1,n)
***********/
  y=0;
  n2=n-n1;
  do i=a+1 to n1;
    y=y+pd(i,p,n1)*(1-cd(c-i,p,n2));
  end;
  return (y);
finish power;

start search(p0,p1) global(alpha,beta);
/*************
Find design (a,c,n1,n) such that
    P(Y1>a,Y>c|p0,n1,n)<=alpha and
    P(Y1>a,Y>c|p1,n1,n)>=1-beta and 
minimize average sample size. The generated design is a Type II, optimal
design with equal sample sizes.

Since we search n from 10 to m, and try to find the satisfactory design with 
smallest average sample size. Thus, the generated design is a optimal design 
(minimizing average sample size).
****************/

  m=100;
  ave=100;
  flag=0;
  do n=10 to m while (ave>n/2);/*search for n*/
	ll=floor(n/2);
	if (2*ll<n) then n1=ll+1;
	else n1=ll;
	n2=n-n1;
    u=a_up(p1,n1);
	do i=0 to u;/*search for a*/
	    aa=(1-cd(i,p0,n1))*n2+n1;
		if (aa<ave) then do;
		  ff=0;
		  cc=i-1;
		  do while ((ff=0) & (cc<n+1));/*search for c*/
		    cc=cc+1;
		    al=power(i,cc,n1,n,p0);
            pw=power(i,cc,n1,n,p1);
		    if ((pw>=1-beta) & (al<=alpha)) then do;
			  ave=aa;
			  flag=1;
			  ff=0;
			  aaa=i;
			  ccc=cc;
			  nnn1=n1;
			  nnn=n;
			  aal=al;
			  ppw=pw;
			  aave=ave;
			  beta_spt=cd(aaa,p1,nnn1);
		    end;
			else if (pw<1-beta) then cc=n+1;
		  end;
	  end;
	end;
  end;
  vv=j(9,1,0);
  vv[1]=flag;
  if flag=1 then do;
    vv[2]=aaa;
    vv[3]=ccc;
    vv[4]=nnn1;
    vv[5]=nnn;
	vv[6]=aal;
	vv[7]=ppw;
	vv[8]=aave;
	vv[9]=beta_spt;
  end;
  return (vv);
finish search;

start print_out(vv);
  a=vv[2];
  c=vv[3];
  n1=vv[4];
  n=vv[5];
  al=vv[6];
  pw=vv[7];
  ave_n=vv[8];
  beta_spt=vv[9];
  print a c n1 n al pw ave_n beta_spt;
finish print_out;

start aa(beta_spt,n1,n,n1_star,p1,ind)global(beta);
/*************************
Given the planned sample sizes of n1 and n, the planned beta spent (beta_spt) at the 
first stage, and the actual sample size of n1_star at the first stage, this subroutine 
finds the largest integer uu such that
            P(Y1<=uu | p1,n1_star)
is closest to btspt, where btspt is the beta spent at the first stage according 
to the sample size n1_star.
ind=1 means the proposed method.
ind=2 means Green's method with fixed spending.
****************************/
if ind=1 then do;
  n2=n-n1;
  if n1_star<=n1 then btspt=beta_spt*n1_star/n1;
  else btspt=beta_spt+(beta-beta_spt)*(n1_star-n1)/n2;
end;
else if ind=2 then do;
  btspt=beta_spt;
end;

  m=-1;
  flag=0;
  do while ((flag=0) & (m<n1_star+1));
    y=cd(m,p1,n1_star);
	if y<=btspt then m=m+1;
	else flag=1;
  end;

  uu=j(2,1,0);
  z=cd(m-1,p1,n1_star);
  if ((y-btspt)>(btspt-z) | (y>beta)) then uu[1]=m-1;
  else uu[1]=m;

  uu[2]=btspt;
  return (uu);
finish aa;

start gen(beta_spt,n1,n,n1_star,n_star,ind) global(alpha,beta,p0,p1);
/***************
For the planned design n1, n, and beta_spt at the first stage, and for the 
attaned sample sizes n1_star and n_star, this subroutine finds the attained design
according to the actual sample sizes n1_star and n_star using type the type II error 
probability spending function.
ind=1 means the proposed method.
ind=2 means Green's method with fixed spending.
***************/ 
  uu=aa(beta_spt,n1,n,n1_star,p1,ind);
  a=uu[1];

  n2_star=n_star-n1_star;
  flag=0;
  c=a-1;
  do while ((flag=0) & (c<n_star+1));
    c=c+1;
	al=power(a,c,n1_star,n_star,p0);
	if al<=alpha then do;
	  flag=1;
	  aaa=a;
	  aal=al;
	  ppw=power(a,c,n1_star,n_star,p1);
	  aave=n2_star*(1-cd(a,p0,n1_star))+n1_star;
	end;
  end;
	vv=j(9,1,0);
	if flag=1 then do;
	  vv[1]=aaa;
	  vv[2]=c;
	  vv[3]=aal;
	  vv[4]=ppw;
	  vv[5]=aave;
      vv[6]=uu[2];
	  vv[7]=cd(aaa,p1,n1_star);
	  vv[8]=n1_star;
	  vv[9]=n_star;
	end;
	return (vv);
  finish gen;

start print1_out(vv);
/***************
Btspt is the computed beta spending at stage 1 using the type II error probability 
spending function defined by the planned design and n1_star. Btspt_star is the actual 
beta spending at stage 1 according to the attained design (a_star, c_star, n1_star, n_star).
**************/
  a_star=vv[1];
  c_star=vv[2];
  al_star=vv[3];
  pw_star=vv[4];
  ave_star=vv[5];
  btspt=vv[6];
  btspt_star=vv[7];
  n1_star=vv[8];
  n_star_star=vv[9];
  print  a_star c_star n1_star n_star_star;
  print al_star pw_star ave_star btspt btspt_star; 
finish print1_out;

start modify(p0,p1,n1_star) global(alpha,beta);
/*************
Find design (a_star,c_star,n1_star,n_star) for given n1_star at stage 1 such that
    P(y1>a_star_star, Y>c_star|p0,n1_star,n_star)<=alpha and
    P(y1>a_star_star, Y>c_star|p1,n1_star,n_star)>=1-beta and
minimizing the average sample size.
****************/
  flag=0;
  ave=100;
  n1=n1_star;
  m=100;
  do n=n1+1 to m;/*search for n*/
	  n2=n-n1;
	  u=a_up(p1,n1);
	  do i=0 to u;/*search for a*/
	    aa=(1-cd(i,p0,n1))*n2+n1;
		if (aa<ave) then do;
		  ff=0;
		  cc=i-1;
	      do while ((ff=0) & (cc<n+1));/*search for c*/
			cc=cc+1;
		    al=power(i,cc,n1,n,p0);
            pw=power(i,cc,n1,n,p1);
            if ((pw>=1-beta) & (al<=alpha)) then do;
			  ave=aa;
		   	  ff=1;
			  flag=1;
			  aaa=i;
			  ccc=cc;
			  nnn1=n1;
			  nnn=n;
			  aal=al;
			  ppw=pw;
			  aave=ave;
			  beta_spt=cd(aaa,p1,nnn1);
		    end;
			else if (pw<1-beta) then cc=n+1;
		  end;
		end;
	  end;
	
  end;
  vv=j(9,1,0);
  vv[1]=flag;
  if flag=1 then do;
    vv[2]=aaa;
    vv[3]=ccc;
    vv[4]=nnn1;
    vv[5]=nnn;
	vv[6]=aal;
	vv[7]=ppw;
	vv[8]=aave;
	vv[9]=beta_spt;
  end;

  a_star=aaa;
  c_star=ccc;
  n_star=nnn;
  pw_star=ppw;
  al_star=aal;
  ave_star=aave;

  print "**********************Modified Design given n1_star**************************";
  print a_star c_star n1_star n_star al_star pw_star ave_star;
  return (vv);
finish modify;

start c_double_star(p0,p1,a_star,n1_star,n_star_star) global(alpha,beta);
/*************
Find c_star_star given a_star, n1_star, and n_star_star such that
    P(y1>a_star, Y>c_star_star|p0,n1_star,n_star_star)<=alpha.
****************/
          n1=n1_star;
		  n=n_star_star;
		  a=a_star;
		  flag=0;
		  ff=0;
		  cc=a-1;
	      do while ((ff=0) & (cc<n+1));/*search for c*/
			cc=cc+1;
		    al=power(a,cc,n1,n,p0);
            pw=power(a,cc,n1,n,p1);
            if (al<=alpha) then do;
			  ave=n1+(1-cd(a,p0,n1_star))*(n-n1);
		   	  ff=1;
			  flag=1;
			  aaa=a;
			  ccc=cc;
			  nnn1=n1;
			  nnn=n;
			  aal=al;
			  ppw=pw;
			  aave=ave;
			  beta_spt=cd(aaa,p1,nnn1);
		    end;
		  end;
		
  vv=j(9,1,0);
  vv[1]=flag;
  if flag=1 then do;
    vv[2]=aaa;
    vv[3]=ccc;
    vv[4]=nnn1;
    vv[5]=nnn;
	vv[6]=aal;
	vv[7]=ppw;
	vv[8]=aave;
	vv[9]=beta_spt;
  end;

  c_star_star=ccc;
  n_star_star=nnn;
  pw_star=ppw;
  al_star=aal;
  ave_star_star=aave;

  print "*********Final Attained Design given n1_star and n_star_star**************************";
  print a_star c_star_star n1_star n_star_star al_star pw_star ave_star_star;

  return (vv);
finish c_double_star;

start main;
  run setup;
  print "Category_I_redesign.out";
  
  do i=1 to 2;
    if i=1 then do;
	  alpha=0.1;
	  beta=0.1;
	end;
    if i=2 then do;
	  alpha=0.05;
	  beta=0.2;
	end;
	do j=1 to 8;
	  if j=1 then do;
	    p0=0.05;
		p1=0.20;
	  end;
	  else do;
	    p0=(j-1)*0.1;
        p1=p0+0.2;
	  end;

      vv=search(p0,p1);
	  Print "////////////////////////////////////New Case/////////////////////////////";
	  print p0 p1 alpha beta;
      run print_out(vv);

	  n1=vv[4];
	  n=vv[5];
	  do k=1 to 2;
	    n1_star=n1+(k-1.5)*4;
		uu=modify(p0,p1,n1_star);
		a_star=uu[2];
		n_star=uu[5];
		do l=1 to 2;
		  n_star_star=n_star+(l-1.5)*4;
          print n1_star n_star n_star_star;
          c_star_star=c_double_star(p0,p1,a_star,n1_star,n_star_star);

          beta_spt=0.02;
          print "**********************Green & Dahlberg***************************";
          print n1_star n_star_star beta_spt;
		  ind=2;
          uu=gen(beta_spt,n1,n,n1_star,n_star_star,ind);
		  run print1_out(uu);
		end;
	  end;
	end;
  end;
finish main;  
 
run;
quit;