/*******************
Category_II_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; if Y1>=b then the drug is claimed as promisssion
and the trial is stopped at the first stage. 
If a<Y1<b 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>=b|p0,n1)+P(a<Y1<b, Y>c|n1,n2,p0).
The power is:              
                     P(y1>=b|p1,n1)+P(a<Y1<b, Y>c|n1,n2,p1).
The average sample size is: 
                     aven=P(a<Y1<b|p0)*n2+n1;
The program search all possible combinaitons of (a,b,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, b_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, b_star, c_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.10;
  beta=0.30;
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|p1,n1)<=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 b_low(p0,n1) global(alpha);
  /***************
Find the smallest integer uu such that
   P(Y1>=uu|p0,n1)<=alpha
***************/
  flag=0;
  m=n1+1;
  do while((flag=0) & (m>0));
    y=1-cd(m-1,p0,n1);
	if y<=alpha then m=m-1;
	else do;
	  flag=1;
	  uu=m+1;
	end;
  end;
  return (uu);
finish b_low;
     
start power(a,b,c,n1,n,p);
/***********
Compute the power of the given design a,b,c,n1,n) under response rate p:
  P(Y1>=b|p,n1)+P(a<y1<b, Y>c|p,n1,n)
***********/
  y=1-cd(b-1,p,n1);
  n2=n-n1;
  do i=a+1 to b-1;
    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,b,c,n1,n) such that
    P(Y1>=b|p0,n1)+P(a<Y1<b,Y>c|p0,n1,n)<=alpha and
    P(Y1>=b|p1,n1)+P(a<Y1<b,Y>c|p1,n1,n)>=1-beta and 
minimize the average sample size. The found 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 
minimum average sample size. Thus, the generated design is a optimal design 
(minimizing average sample size).
****************/
  m=100;
  ave=100;
  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);
	v=b_low(p0,n1);
	do i=0 to u;/*search for a*/
	  do j=v to n1;/*search for b*/
	    aa=(cd(j-1,p0,n1)-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,j,cc,n1,n,p0);
            pw=power(i,j,cc,n1,n,p1);
		    if ((pw>=1-beta) & (al<=alpha)) then do;
			  ave=aa;
			  flag=1;
			  ff=0;
			  aaa=i;
			  bbb=j;
			  ccc=cc;
			  nnn1=n1;
			  nnn=n;
			  aal=al;
			  ppw=pw;
			  aave=ave;
			  beta_spt=cd(aaa,p1,nnn1);
			  alpha_spt=1-cd(bbb-1,p0,nnn1);
		    end;
			else if (pw<1-beta) then cc=n+1;
		  end;
		end;
	  end;
	end;
  end;
  vv=j(11,1,0);
  vv[1]=flag;
  if flag=1 then do;
    vv[2]=aaa;
	vv[3]=bbb;
    vv[4]=ccc;
    vv[5]=nnn1;
    vv[6]=nnn;
	vv[7]=aal;
	vv[8]=ppw;
	vv[9]=aave;
	vv[10]=beta_spt;
	vv[11]=alpha_spt;
  end;
  return (vv);
finish search;

start print_out(vv);
  a=vv[2];
  b=vv[3];
  c=vv[4];
  n1=vv[5];
  n=vv[6];
  al=vv[7];
  pw=vv[8];
  ave_n=vv[9];
  beta_spt=vv[10];
  alpha_spt=vv[11];
  print a b c n1 n al pw ave_n beta_spt alpha_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 bb(alpha_spt,n1,n,n1_star,p0,ind) global(alpha);
/*************************
Given the planned sample sizes of n1 and n, the planned alpha spent (alpha_spt) at the 
first stage, and the actual sample size of n1_star at the first stage, this subroutine 
finds the smallest integer uu such that
            P(Y1>=uu | p1,n1_star)
is closest to alspt, where alspt is the alpha 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 alspt=alpha_spt*n1_star/n1;
  else alspt=alpha_spt+(alpha-alpha_spt)*(n1_star-n1)/n2;
end;
else if ind=2 then do;
  alspt=alpha_spt;
end;

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

  uu=j(2,1,0);
  z=1-cd(m,p0,n1_star);
  if ((y-alspt)>(alspt-z) | (y>alpha)) then uu[1]=m+1;
  else uu[1]=m;

  uu[2]=alspt;
  return (uu);
finish bb;

start gen(beta_spt,alpha_spt,n1,n,n1_star,n_star,ind) global(alpha,beta,p0,p1);
/***************
For the planned design n1, n, beta_spt, and alpha_spt at the first stage, and for the 
attaned sample sizes n1_star and n_star, this subroutine finds the design
according to the actual sample sizes n1_star and n_star using type I and type II error 
probability spending functions.
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);
  ww=bb(alpha_spt,n1,n,n1_star,p0,ind);
  a=uu[1];
  b=ww[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,b,c,n1_star,n_star,p0);
	if al<=alpha then do;
	  flag=1;
	  aaa=a;
	  bbb=b;
	  aal=al;
	  ppw=power(a,b,c,n1_star,n_star,p1);
	  aave=n2_star*(cd(b-1,p0,n1_star)-cd(a,p0,n1_star))+n1_star;
	end;
  end;
	vv=j(12,1,0);
	if flag=1 then do;
	  vv[1]=aaa;
	  vv[2]=bbb;
	  vv[3]=c;
	  vv[4]=aal;
	  vv[5]=ppw;
	  vv[6]=aave;
      vv[7]=uu[2];
      vv[8]=ww[2];
	  vv[9]=cd(aaa,p1,n1_star);
	  vv[10]=1-cd(bbb-1,p0,n1_star);
	  vv[11]=n1_star;
	  vv[12]=n_star;
	end;
	return (vv);
  finish gen;

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

start modify(p0,p1,n1_star) global(alpha,beta);
/*************
Find design (a_star,b_star,c_star,n1_star,n_star) for given n1_star such that
    P(Y1>=b_star|p0,n1_star)+P(a_star<y1<b_star, Y>c_star|p0,n1_star,n_star)<=alpha and
    P(Y1>=b_star|p0,n1_star)+P(a_star<y1<b_star, Y>c_star|p0,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);
	  v=b_low(p0,n1);
	  do i=0 to u;/*search for a*/
	  do j=v to n1;/*search for b*/
	    aa=(cd(j-1,p0,n1)-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,j,cc,n1,n,p0);
            pw=power(i,j,cc,n1,n,p1);
            if ((pw>=1-beta) & (al<=alpha)) then do;
			  ave=aa;
		   	  ff=1;
			  flag=1;
			  aaa=i;
			  bbb=j;
			  ccc=cc;
			  nnn1=n1;
			  nnn=n;
			  aal=al;
			  ppw=pw;
			  aave=ave;
			  beta_spt=cd(aaa,p1,nnn1);
			  alpha_spt=1-cd(bbb-1,p0,nnn1);
		    end;
			else if (pw<1-beta) then cc=n+1;
		  end;
		end;
	  end;
	  end;
	
  end;
  vv=j(11,1,0);
  vv[1]=flag;
  if flag=1 then do;
    vv[2]=aaa;
	vv[3]=bbb;
    vv[4]=ccc;
    vv[5]=nnn1;
    vv[6]=nnn;
	vv[7]=aal;
	vv[8]=ppw;
	vv[9]=aave;
	vv[10]=beta_spt;
	vv[11]=alpha_spt;
  end;

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

  print "**********************Modified Design given n1_star**************************";
  print a_star b_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,b_star,n1_star,n_star_star) global(alpha,beta);
/*************
Find c_star_star given a_star, b_star n1_star, and n_star_star such that
    P(Y1>=b_star|p0,n1_star)+P(a_star<y1<b_star, Y>c_star_star|p0,n1_star,n_star_star)<=alpha.
****************/
          n1=n1_star;
		  n=n_star_star;
		  a=a_star;
		  b=b_star;
		  flag=0;
		  ff=0;
		  cc=a-2;
	      do while ((ff=0) & (cc<n+1));/*search for c*/
			cc=cc+1;
		    al=power(a,b,cc,n1,n,p0);
            pw=power(a,b,cc,n1,n,p1);
            if (al<=alpha) then do;
			  ave=n1+(cd(b-1,p0,n1)-cd(a,p0,n1_star))*(n-n1);
		   	  ff=1;
			  flag=1;
			  aaa=a;
			  bbb=b;
			  ccc=cc;
			  nnn1=n1;
			  nnn=n;
			  aal=al;
			  ppw=pw;
			  aave=ave;
			  beta_spt=cd(aaa,p1,nnn1);
			  alpha_spt=1-cd(bbb-1,p0,nnn1);
		    end;
		  end;
		
  vv=j(11,1,0);
  vv[1]=flag;
  if flag=1 then do;
    vv[2]=aaa;
	vv[3]=bbb;
    vv[4]=ccc;
    vv[5]=nnn1;
    vv[6]=nnn;
	vv[7]=aal;
	vv[8]=ppw;
	vv[9]=aave;
	vv[10]=beta_spt;
	vv[11]=alpha_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 b_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_II_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[5];
	  n=vv[6];
	  do k=1 to 2;
	    n1_star=n1+(k-1.5)*4;
		uu=modify(p0,p1,n1_star);
		a_star=uu[2];
		b_star=uu[3];
		n_star=uu[6];
		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,b_star,n1_star,n_star_star);

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