* Design of Group Sequential Trials for Simultaneous Efficacy and Futility Testing;

%macro EffFutDesign(fraction, effsize, power, alpha, rhoeff, rhofut, boundary, sizepower);
/*
Inputs:
 
FRACTION = Input data set that contains fractions of the total sample 
           size accrued at successive analyses

EFFSIZE  = True effect size
 
POWER    = Power

ALPHA    = One-sided Type I error probability

RHOEFF   = Shape parameter of upper (efficacy) stopping boundary
           (0.5 if Pocock boundary and 0 if O'Brien-Fleming boundary)

RHOFUT   = Shape parameter of lower (futility) stopping boundary
           (0.5 if Pocock boundary and 0 if O'Brien-Fleming boundary)

BOUNDARY = Output data set that contains stopping probabilities at scheduled 
           looks

SIZEPOWER= Output data set that contains average sample number and power for 
           selected effect sizes

*/
proc iml;
    start ParSearch(c) global(m,lastlook,stfract,inc);
        drift=(c[1]*lastlook##&rhoeff+c[2]*lastlook##&rhofut)/lastlook;
        upper=c[1]*stfract##&rhoeff;
        lower=drift*stfract-c[2]*stfract##&rhofut;
        lower[m]=lower[m]-1e-5;
        boundary=lower//upper;
        call seq(p,boundary) eps=1e-8 tscale=inc;
        crossh0=sum(p[3,]-p[2,])-α
        adjustment=drift*stfract;
        boundary=(lower-adjustment)//(upper-adjustment);
        call seq(p,boundary) eps=1e-8 tscale=inc;
        crossh1=sum(p[3,]-p[2,])-&power;
        diff=abs(crossh0)+abs(crossh1);
        return(diff);
    finish;
    use &fraction;
    read all var _all_ into fraction;    
    m=nrow(fraction);    
    fract=t(fraction);    
    stfract=fract/fract[1];   
    inc=j(1,m-1,0);
    do i=1 to m-1;
        inc[i]=(fract[i+1]-fract[i])/fract[1];
    end;     
    lastlook=stfract[m];
    nfixed=2*((probit(&power)+probit(1-&alpha))/&effsize)**2;      
    start={1 1};
    tc=repeat(.,1,12);
    tc[1]=100;
    tc[3]=1e-5;
    call nlpdd(rc,c,"ParSearch",start) tc=tc;
    drift=(c[1]*lastlook##&rhoeff+c[2]*lastlook##&rhofut)/lastlook;
    upper=c[1]*stfract##&rhoeff;
    lower=drift*stfract-c[2]*stfract##&rhofut;    
    lower[m]=lower[m]-1e-5;
    boundary=lower//upper;
    call seq(prob0,boundary) eps=1e-8 tscale=inc;    
    adjustment=drift*stfract;
    boundary=(lower-adjustment)//(upper-adjustment);               
    call seq(prob1,boundary) eps=1e-8 tscale=inc;
    upperz=(stfract##(-0.5))#upper;
    lowerz=(stfract##(-0.5))#lower;        
    upperp=1-probnorm(upperz);
    lowerp=1-probnorm(lowerz);
    max=2*(drift/&effsize)*(drift/&effsize)/fract[1];        
    boundary=j(m,10,0);
    boundary[,1]=t(1:m);
    boundary[,2]=ceil(fraction*max);
    boundary[,3]=t(lowerz);
    boundary[,4]=t(upperz);
    boundary[,5]=t(lowerp);
    boundary[,6]=t(upperp);
    boundary[,7]=t(prob0[3,]-prob0[2,]+prob0[1,]); 
    boundary[,8]=t(cusum(prob0[3,]-prob0[2,]+prob0[1,]));     
    boundary[,9]=t(prob1[3,]-prob1[2,]+prob1[1,]); 
    boundary[,10]=t(cusum(prob1[3,]-prob1[2,]+prob1[1,]));             
    varnames={"Analysis", "Size", "LowerTestStBoundary", "UpperTestStBoundary",
        "LowerPValBoundary", "UpperPValBoundary", 
        "ProbH0", "CumProbH0", "ProbH1", "CumProbH1"};    
    create &boundary from boundary[colname=varnames];
    append from boundary;  
    sizepower=j(21,3,0);
    do i=0 to 20;          
        adjustment=i*drift*stfract/10;
        boundary=(lower-adjustment)//(upper-adjustment);           
        call seq(prob2,boundary) eps=1e-8 tscale=inc; 
        stop=prob2[3,]-prob2[2,]+prob2[1,]; 
        sizepower[i+1,1]=i*&effsize/10;
        sizepower[i+1,2]=ceil(max*(1-(1-fract)*stop`));   
        sizepower[i+1,3]=sum(prob2[3,]-prob2[2,]);*1-(prob2[2,]-prob2[1,])[m]; 
    end;
    varnames={"EffSize", "AveSize", "Power"};    
    create &sizepower from sizepower[colname=varnames];
    append from sizepower;        
    summary=j(1,4,0);
    summary[1]=ceil(max); summary[2]=sizepower[1,2]; summary[3]=sizepower[11,2];    
    summary[4]=ceil(nfixed);
    create summary from summary;
    append from summary;        
    quit;
data summary;
    set summary;
    format value best6.;
    length par $50;
    par="One-sided Type I error probability"; value=α output;    
    par="Power"; value=&power; output;
    par="True effect size"; value=&effsize; output;
    par="Shape parameter of upper boundary"; value=&rhoeff; output;
    par="Shape parameter of lower boundary"; value=&rhofut; output;
    par="Maximum sample size per group"; value=col1; output;
    par="Average sample size per group under H0"; value=col2; output;
    par="Average sample size per group under H1"; value=col3; output;
    par="Fixed sample size per group"; value=col4; output;
    label par="Summary" value="Value";
    keep par value;
proc print data=summary noobs label;
    var par value;
data &boundary;
    set &boundary;
    format LowerTestStBoundary UpperTestStBoundary LowerPValBoundary 
        UpperPValBoundary ProbH0 CumProbH0 ProbH1 CumProbH1 6.4 
        Analysis Size 4.0;
    label Analysis="Analysis"
          Size="Sample size per group"
          LowerTestStBoundary="Lower stopping boundary (test statistic scale)"
          UpperTestStBoundary="Upper stopping boundary (test statistic scale)"
          LowerPValBoundary="Lower stopping boundary (p-value scale)"
          UpperPValBoundary="Upper stopping boundary (p-value scale)"
          ProbH0="Stopping probability under H0"
          CumProbH0="Cumulative stopping probability under H0"
          ProbH1="Stopping probability under H1"
          CumProbH1="Cumulative stopping probability under H1";
data &sizepower;
    set &sizepower;
    format EffSize best6. AveSize 5.0;
    label EffSize="True effect size"
          AveSize="Average sample size per group"
          Power="Power";
    run;    
%mend EffFutDesign;