/*******************
  Improved two-stage tests.sas.

  This program generates designs as in Table I in the paper entitled 
  "Improved two-stage tests for stratified phase II cancer clinical trials",
  Statistics in Medicine, 2012, by Chang, MN, JJ Shuster, and W hou.
 
 
  The number of strata is k = 3. The input parameters are p01, p11, p02, p12, 
  delta1, delta2, delta3, alpha, beta (1 - power), afract (gamma1), and bfrat (gamma2)
  are described in the paper, where p11 = p10 + delta1, p21=p20 + delta2, and 
  p31= p30 +delta 3. 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 I 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;
  
  p10=0.5;
  p20=0.3;
  p30=0.1;
  delta1=0.2;
  delta2=0.2;
  delta3=0.2;
  nn11=2;
  nn12=4;
  nn13=14;
  alpha=0.05;
  beta=0.2;
  afract=0.25;
  bfract=0.25;

  epsilon=10**(-8);
  zz=j(3,1);
  ww=j(3,1);
  zz[1]=p10;  ww[1]=delta1;
  zz[2]=p20;  ww[2]=delta2;
  zz[3]=p30;  ww[3]=delta3;
finish setup;

start rej2(n1,n2,p1,p2,w1,w2,b) global(epsilon);
/****
Compute PR(w1*Y1+w2*Y2>=b), where Y1 and Y2 ar independent binomial
random variables with parameters(n1,p1) and (n2,p2);
****/
 * print n1 n2 p1 p2 w1 w2 b;
  y=0;
  b=b-epsilon;
  if b<=0 then y=1;
  else if b>w1*n1+w2*n2 then y=0;
  else if n2=0 then do;
    x=int(b/w1);
    y=1-cdf('binomial',x,p1,n1);
  end;
  else do;
    y=0;
    do i=0 to n1;
	  x=int((b-w1*i)/w2);
  *print i p1 n1 b w1 w2 p2 n2 x;
	  if (b-w1*i)<0 then
	    y=y+pdf('binomial',i,p1,n1);
	  else
        y=y+pdf('binomial',i,p1,n1)*(1-cdf('binomial',x,p2,n2));
    end;
  end;
  return (y);
finish rej2;

start rej3(n1,n2,n3,p1,p2,p3,w1,w2,w3,b) global(epsilon);
/****
Compute PR(w1*Y1+w2*Y2+w3*Y3>=b), where Y1, Y2, and Y3 are independent binomial
random variables with parameters(n1,p1),(n2,p2), and (n3,p3);
****/
  y=0;
  if b<=0 then y=1;
  else if b>w1*n1+w2*n2+w3*n3 then y=0;
  else if n3=0 then
    y=rej2(n1,n2,p1,p2,w1,w2,b);
  else do;
    y=0;
    do i=0 to n1;
      y=y+pdf('binomial',i,p1,n1)*rej2(n2,n3,p2,p3,w2,w3,b-w1*i);
    end;
  end;
  return (y);
finish rej3;

start bdeterm(t) global(nn11,nn12,nn13,w1,w2,w3,zz,alpha,afract,epsilon);
/*****************
Determine the upper bound b at the first stage. If R1 = w1*y11+w2*y12+w3*y13 > b 
then reject H0.
******************/
  a=alpha*afract;
  p1=zz[1]; p2=zz[2]; p3=zz[3]; 
  y=1;
  aa=-1;
  bb=int(w1*nn11+w2*nn12+w3*nn13)+1;
  do until (bb-aa<0.00001);
    c=(aa+bb)/2;
    y=rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,c+epsilon);
	if y>a then aa=c;
	else bb=c;
  end;
/******
  u=rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,aa+epsilon)-a;
  v=a-rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,bb+epsilon);
  if u<v then w=aa;
  else w=bb;
********/
  bb=bb+0.001;
  return (bb);
finish bdeterm;

start adeterm(t) global(nn11,nn12,nn13,w1,w2,w3,zz,ww,beta,bfract,epsilon);
/*****************
Determine the lower bound a at the first stage. If R1 = w1*y11+w2*y12+w3*y13 <a 
then accept H0.
******************/
  b=beta*bfract;
  p1=zz[1]+ww[1]; p2=zz[2]+ww[2]; p3=zz[3]+ww[3]; 
  y=0;
  aa=-1;
  bb=int(w1*nn11+w2*nn12+w3*nn13)+1;
  do until (bb-aa<0.00001);
    c=(aa+bb)/2;
    y=1-rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,c);
*print aa bb c y;
	if y<=b then aa=c;
	else bb=c;
*print aa bb c y;
  end;
/*********
  u=b-1+rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,aa);
  v=1-rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,bb)-b;
  if u<v then w=aa;
  else w=bb;

u=1-rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,aa);
v=1-rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,bb);
print aa bb u v;
*************/
  aa=aa-0.001;
  return (aa);
finish adeterm;

start alp(a,b,c) global(nn11,nn12,nn13,nn21,nn22,nn23,w1,w2,w3,zz,epsilon);
/********************
With the upper bound b and the lower bound a at the first stage 
and the bound c at the second stage, compute the significance level 
of the test: reject H0 if R1 = w1*y11+w2*y12+w3*y13 >b; accept H0 if R1<a; 
continue to the second stage when a<=R1<=b and reject H0 if R1+R2>c.
**********************/
  p1=zz[1]; p2=zz[2]; p3=zz[3]; 
  c=c+epsilon;
  z=rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,b);
  do i1=0 to nn11;
    do i2=0 to nn12;
	  do i3=0 to nn13;
	    if (a<=w1*i1+w2*i2+w3*i3) & (w1*i1+w2*i2+w3*i3<=b) then do;
          y=pdf('binomial',i1,p1,nn11)*pdf('binomial',i2,p2,nn12)*
            pdf('binomial',i3,p3,nn13);
          u=rej3(nn21,nn22,nn23,p1,p2,p3,w1,w2,w3,c-w1*i1-w2*i2-w3*i3);
          z=z+y*u;
        end;
	  end;
    end;
  end;
  return (z);
finish alp;

start bet(a,b,c) global(nn11,nn12,nn13,nn21,nn22,nn23,w1,w2,w3,zz,ww,epsilon);
/********************
With the upper bound b and the lower bound a at the first stage 
and the bound c at the second stage, compute the power
of the test: reject H0 if R1 = w1*y11+w2*y12+w3*y13 >b; accept H0 if r1<a; continue to 
the second stage when a<=R1<=b and reject H0 if R1+R2>c.
**********************/
  c=c+epsilon;
  p1=zz[1]+ww[1]; p2=zz[2]+ww[2]; p3=zz[3]+ww[3]; 
  z=1-rej3(nn11,nn12,nn13,p1,p2,p3,w1,w2,w3,a);
  do i1=0 to nn11;
    do i2=0 to nn12;
	  do i3=0 to nn13;
	    if (a<=w1*i1+w2*i2+w3*i3) & (w1*i1+w2*i2+w3*i3<=b) then do;
          y=pdf('binomial',i1,p1,nn11)*pdf('binomial',i2,p2,nn12)*
            pdf('binomial',i3,p3,nn13);
          u=1-rej3(nn21,nn22,nn23,p1,p2,p3,w1,w2,w3,c-w1*i1-w2*i2-w3*i3);
          z=z+y*u;
         end;
       end;
	 end;
  end;
  return (z);
finish bet;

start cdeterm(a,b) global(nn11,nn12,nn13,nn21,nn22,nn23,w1,w2,w3,zz,alpha,epsilon);
/**********************
Giben the boundaries at the first stage ( a and b), determine the bound 
at the second stage to satisfy the significance requirement.
***********************/
  p1=zz[1]; p2=zz[2]; p3=zz[3]; 
  aa=-1;
  bb=int(w1*(nn11+nn21)+w2*(nn12+nn22)+w3*(nn13+nn23))+1;
  do until (bb-aa<0.00001);
    c=(aa+bb)/2;
	y=alp(a,b,c);
	if y<=alpha then bb=c;
	else aa=c;
  end;
  bb=bb+0.001;
  return (bb);
finish cdeterm;

start prt;

  print nn11 nn12 nn13;
  print alpha beta delta1 delta2 delta3;
  print p10 p20 p30;

  b=bdeterm(1);
  a=adeterm(1);
  c=cdeterm(a,b);
  a=(floor(a*1000)+1)/1000;
  b=(floor(b*1000))/1000;
  c=(floor(c*1000))/1000;
  al=alp(a,b,c);
  bt=bet(a,b,c);
  pw=1-bt;
  print a b c al pw;

finish prt;

start main;
  run setup;

    zz[1]=p10;  ww[1]=delta1; 
    zz[2]=p20;  ww[2]=delta2;
    zz[3]=p30;  ww[3]=delta3;
	p11=p10+delta1;
    p21=p20+delta2;
    p31=p30+delta3;

    w1=log(p11*(1-p10)/(1-p11)/p10);
	w2=log(p21*(1-p20)/(1-p21)/p20);
    w3=log(p31*(1-p30)/(1-p31)/p30);

    nn21=nn11;nn22=nn12;nn23=nn13;

    run prt;
   
finish main;

run;
quit;