* Design of Group Sequential Trials for Simultaneous Efficacy Testing;



%macro EffDesign(fraction, effsize, power, alpha, rho, 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

RHO      = Shape parameter of 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 DriftSearch(d) global(m,critvalue,stfract,inc,infinity);        
        upper=critvalue*stfract##ρ
        adjustment=d*stfract;
        boundary=infinity//(upper-adjustment);                     
        call seq(p,boundary) eps=1e-8 tscale=inc; 
        diff=abs(1-&power-(p[2,]-p[1,])[m]);          
        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;        
    infinity=repeat(.m,1,m);
    upper=stfract##ρ
    boundary=infinity//upper;    
    call seqscale(prob,critvalue,boundary,1-&alpha) eps=1e-8 tscale=inc; 
    upper=critvalue*stfract##ρ
    boundary=infinity//upper;    
    stopz=critvalue*stfract##(&rho-0.5);    
    stopp=1-probnorm(stopz);
    call seq(prob0,boundary) eps=1e-8 tscale=inc; 
    nfixed=2*((probit(&power)+probit(1-&alpha))/&effsize)**2;    
    start=&effsize*sqrt(nfixed*fract[1]/2);
    tc=repeat(.,1,12);
    tc[1]=100;
    tc[3]=1e-5;
    call nlpdd(rc,drift,"DriftSearch",start) tc=tc;
    max=2*(drift/&effsize)*(drift/&effsize)/fract[1];
    upper=critvalue*stfract##ρ    
    adjustment=drift*stfract;
    boundary=infinity//(upper-adjustment);                 
    call seq(prob1,boundary) eps=1e-8 tscale=inc;   
    &boundary=j(m,8,0);
    &boundary[,1]=t(1:m);
    &boundary[,2]=ceil(fraction*max);
    &boundary[,3]=t(stopz);
    &boundary[,4]=t(stopp);
    &boundary[,5]=t(prob0[3,]-prob0[2,]+prob0[1,]); 
    &boundary[,6]=t(cusum(prob0[3,]-prob0[2,]+prob0[1,]));     
    &boundary[,7]=t(prob1[3,]-prob1[2,]+prob1[1,]); 
    &boundary[,8]=t(cusum(prob1[3,]-prob1[2,]+prob1[1,]));             
    varnames={"Analysis", "Size", "TestStBoundary", "PValBoundary", "ProbH0", 
        "CumProbH0", "ProbH1", "CumProbH1"};    
    create &boundary from &boundary[colname=varnames];
    append from &boundary;    
    &sizepower=j(21,3,0);
    do i=0 to 20;
        upper=critvalue*stfract##ρ
        adjustment=i*&effsize*sqrt(max*fract[1]/2)*stfract/10;
        boundary=infinity//(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]=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="Stopping boundary parameter"; value=ρ 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 TestStBoundary PValBoundary ProbH0 CumProbH0 ProbH1 CumProbH1 6.4
        Analysis Size 4.0;
    label Analysis="Analysis"
          Size="Sample size per group"
          TestStBoundary="Stopping boundary (test statistic scale)"
          PValBoundary="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 EffDesign;