/************************************************************************   */
/* These Computer Programs are copyright protected and can only             */
/* be used or modified by those who own the book:                           */
/* "Adaptive Design Theory and Implementation Using SAS and R."             */
/* When you use or make modifications, please do not remove          */
/* the reference to the book. Any comments will be appreciated.             */
/* Thank you for using them.                                                */
/*                     -- Mark.Chang@statisticians.org                       */
/************************************************************************   */
/************************************************************************   */
 

Corrections and improvements for the computer programs (in Bold) have been made since the first print of the book. Some of them were discovered by readers. Thanks Dr. Tad Archambault and Tracy Zhang for pointing out the errors in SAS Macros 2.4 and 6.3. Other changes were made in SAS macros 15.1 and 15.2, and R-function B.6. All those changes do not affect the results of the examples presented in the book, except the Table 15.2, which is corrected in the 2nd print of the book. An improvement has been made in SAS macro 9.1, which leads a higher target conditional power (>95%) for the results in Tables 9.2-9.5.

 
/*
List of SAS Macros
    SAS Macro 2.1: Equivalence Trial with Normal Endpoint
    SAS Macro 2.2: Equivalence Trial with Binary Endpoint
    SAS Macro 2.3: Crossover Bioequivalence Trial
    SAS Macro 2.4: Sample-Size for Dose-Response Trial
    SAS Macro 4.1: Two-Stage Adaptive Design with Binary Endpoint
    SAS Macro 4.2: Two-Stage Adaptive Design with Normal Endpoint
    SAS Macro 4.3: Two-Stage Adaptive Design with Survival Endpoint
    SAS Macro 4.4: Event-Based Adaptive Design
    SAS Macro 4.5: Adaptive Equivalence Trial Design
    SAS Macro 5.1: Stopping Boundaries with Adaptive Designs
    SAS Macro 5.2: Two-Stage Design with Inverse-Normal Method
    SAS Macro 6.1: N-Stage Adaptive Designs with Normal Endpoint
    SAS Macro 6.2: N-Stage Adaptive Designs with Binary Endpoint
    SAS Macro 6.3: N-Stage Adaptive Designs with Various Endpoint
    SAS Macro 7.1: Conditional Power
    SAS Macro 7.2: Sample-Size Based on Conditional Power
    SAS Macro 9.1: Adaptive Design with Sample-Size Reestimation
    SAS Macro 11.1: Two-Stage Drop-Loser Adaptive Design
    SAS Macro 12.1: Biomarker-Adaptive Design
    SAS Macro 14.1: Randomized Play-the-Winner Design
    SAS Macro 14.2: Binary Response-Adaptive Randomization
    SAS Macro 14.3: Normal Response-Adaptive Randomization
    SAS Macro 15.1: 3 + 3 Dose-Escalation Design
    SAS Macro 15.2: Continual Reassessment Method
    SAS Macro 16.1: Simon Two-Stage Futility Design
    SAS Macro A.1: Mixed Exponential Distribution
    SAS Macro A.2: Multi-Variate Normal Distribution
    
               ======================================================*/
 
 
        /*››SAS Macro 2.1: Equivalence Trial with Normal Endpoint››*/
 
    %Macro EquivCI(nSims=1000, nPerGrp=200, ux=0, uy=1, delta=1.2, 
          sigmax=1, sigmay=1.2, alpha=0.05);
    Data TwoGVars;
    Keep xMean yMean powerCI powerTest;
    powerCI=0; powerTest=0;
    Do iSim=1 To &nSims;
     xMean=0; yMean=0; s2x=0; s2y=0;
     Do iObs=1 To &nPerGrp;
        xNOR=Rannor(7362); xMean=xMean+xNor; s2x=s2x+xNor**2;
        yNOR=Rannor(2637); yMean=yMean+yNor; s2y=s2y+yNor**2;
     End;
     xMean=xMean*&sigmax/&nPerGrp+&ux;
     yMean=yMean*&sigmay/&nPerGrp+&uy;
     sp=((s2x*&sigmax**2+s2y*&sigmay**2)/(2*&nPerGrp-2))**0.5;
     se=sp/(&nPerGrp/2)**0.5;
     * CI method;
     ICW=Probit(1-&alpha)*se;
     If Abs(yMean-xMean)+ICW < &delta Then 
          powerCI=powerCI+1/&nSims;
          *Two one-sided test method;
     T1=(xMean-yMean-&delta)/se;
     T2=(xMean-yMean+&delta)/se;
     If T1<-Probit(1-&alpha) & T2>Probit(1-&alpha) Then 
          powerTest=powerTest+1/&nSims; 
    End;
    Output;
    Run;
    Proc Print Data=TwoGVars(obs=1); Run;
    %Mend EquivCI;
  
    
    /*  SAS Macro 2.2: Equivalence Trial with Binary Endpoint*/
 
    %Macro TwoSamZTest(nSims=100000, nPerGrp=100,
             px=0.3, py=0.4, delta=0.3, alpha=0.05);
    Data TwoGVars;
    KEEP powerCI powerTest;
    powerCI=0; powerTest=0;
    Do iSim=1 To &nSims;
     PxObs=Ranbin(733,&nPerGrp,&px)/&nPerGrp;          
     PyObs=Ranbin(236,&nPerGrp,&py)/&nPerGrp;             
     se=((PxObs*(1-PxObs)+PyObs*(1-PyObs))/&nPerGrp)**0.5;
     *CI method;
     ICW=Probit(1-&alpha)*se;
     IF Abs(PxObs-PyObs)+ICW < &delta Then 
             powerCI=powerCI+1/&nSims;
     *Two one-sided test method;
     T1=(PyObs-PxObs-&delta)/se;
     T2=(PyObs-PxObs+&delta)/se;
     IF T1<-Probit(1-&alpha) & T2>Probit(1-&alpha) Then 
             powerTest=powerTest+1/&nSims; 
    End;
    Output;
    Run;
    Proc Print; Run;
    %Mend TwoSamZTest;
 
        /*SAS Macro 2.3: Crossover Bioequivalence Trial*/
 
    %Macro Power2By2ABE(totalN=24, sWithin=0.355, uRatio=1);
    Data ABE; Keep sWithin uRatio n power;
    n=&totalN; sWithin=&sWithin; uRatio=&uRatio;
    * Err df for AB/BA crossover design;
    n2=n-2;
    t1=tinv(1-0.05,n-2); t2=-t1;
    nc1=Sqrt(n)*log(uRatio/0.8)/Sqrt(2)/sWithin;
    nc2=Sqrt(n)*log(uRatio/1.25)/Sqrt(2)/sWithin;
    df=Sqrt(n-2)*(nc1-nc2)/(2*t1);
    Power=Probt(t2,df,nc2)-Probt(t1,df,nc1);
    Run;
    Proc Print; Run;
    %Mend Power2By2ABE;
 
        /*SAS Macro 2.4: Sample-Size for Dose-Response Trial*/
 
    %Macro AsympN(endpoint="normal", alpha=0.025, power=0.8,
                 nArms=5, delta=0, tStd=12, tAcr=4);
    Data contrastN; Set dInput; 
    Keep Endpoint nArms alpha power TotalSampleSize;
    Array u{&nArms}; Array s{&nArms}; Array f{&nArms}; 
    Array c{&nArms}; endpoint=&endpoint; delta=&delta; 
    alpha=&alpha; power=&power; nArms=&nArms;
     epi = 0; s0 = 0;
     Do i =1 To nArms; epi = epi + c{i}*u{i}- &delta/nArms; End;
     If &endpoint = "normal" Then Do;
       Do i =1 To nArms; s0 = s0 + s{i}/nArms; End;
     End;
     If &endpoint = "binary" Then Do;
       Do i = 1 To nArms; 
          s{i} = (u{i}*(1-u{i}))**0.5;
          s0=s0 + s{i}/nArms;
       End;
     End;
     If &endpoint = "survival" Then Do;
        Do i = 1 To nArms;
          s{i} = u{i}*(1+exp(-u{i}*&tStd)*(1-exp(u{i}*&tAcr))
             /(&tAcr*u{i}))**(-0.5);
          s0 = s0 + s{i}/nArms;
        End;
     End;
     sumscf0 = 0; sumscf = 0;
     Do i = 1 To nArms;  sumscf0 = sumscf0 + s0**2*c{i}*c{i}/f{i}; End;
     Do i = 1 To nArms;  sumscf = sumscf + s{i}**2*c{i}*c{i}/f{i}; End;
     n = ((PROBit(1-&alpha)*sumscf0**0.5
          + Probit(&power)*sumscf**0.5)/epi)**2;
     TotalSampleSize = round(n);
    run;
    proc print;
    run;
    %Mend AsympN;
 
        /*SAS Macro 4.1: Two-Stage Adaptive Design with Binary Endpoint*/
 
    %Macro DCSPbinary(nSims=1000000, Model="sum", alpha=0.025,
              beta=0.2, NId=0, Px=0.2, Py=0.4, DuHa=0.2, 
             nAdj="N", Nmax=100, N0=100, nInterim=50, a=2, 
             alpha1=0.01, beta1=0.15, alpha2=0.1871);
    Data DCSPbinary; Keep Model FSP ESP AveN Power nClassic;
      seedx=2534; seedy=6762; Model=&model; NId=&NId; 
     Nmax=&Nmax; N1=&nInterim; Px=&Px; Py=&Py; 
      eSize=abs((&DuHa+NId)/((Px*(1-Px)+Py*(1-Py))/2)**0.5);
      nClassic=Round(2*((probit(1-&alpha)+probit(1-&beta))/eSize)**2);
      FSP=0; ESP=0; AveN=0; Power=0;
      Do isim=1 To &nSims;
      nFinal=N1;
       Px1=Ranbin(seedx,N1,px)/N1;
       Py1=Ranbin(seedy,N1,py)/N1;
       sigma=((Px1*(1-Px1)+Py1*(1-Py1))/2)**0.5; 
       T1 = (Py1-Px1+NId)*Sqrt(N1/2)/sigma;
       p1=1-ProbNorm(T1); 
       If p1>&beta1 Then FSP=FSP+1/&nSims;
       If p1<=&alpha1 Then Do;
          Power=Power+1/&nSims; ESP=ESP+1/&nSims;
       End;
       If p1>&alpha1 and p1<=&beta1 Then Do;
          eRatio=abs(&DuHa/(abs(Py1-Px1)+0.0000001));
          nFinal=Min(&Nmax,Max(&N0,eRatio**&a*&N0)); 
          If &nAdj="N" then nFinal=&Nmax;
          If nFinal>N1 Then Do;
             N2=nFinal-N1;
             Px2=Ranbin(seedx,N2,px)/N2;
             Py2=Ranbin(seedy,N2,py)/N2;
             sigma=((Px2*(1-Px2)+Py2*(1-Py2))/2)**0.5;
             T2 = (Py2-Px2+NId)*Sqrt(N2/2)/sigma;
             p2=1-ProbNorm(T2); 
             If Model="ind" Then TS2=p2; 
             If Model="sum" Then TS2=p1+p2;
             If Model="prd" Then TS2=p1*p2;
             If .<TS2<=&alpha2 then Power=Power+1/&nSims;
          End;
       End;
       AveN=Aven+nFinal/&nSims;
     End;
     Output;
    Run;
    Proc Print Data=DCSPbinary; Run; 
    %Mend DCSPbinary;
 
        /*SAS Macro 4.2: Two-Stage Adaptive Design with Normal Endpoint*/
 
    %Macro DCSPnormal(nSims=1000000, Model="sum", alpha=0.025, 
          beta=0.2, sigma=2, NId=0, ux=0, uy=1, nInterim=50, 
          Nmax=100, N0=100, DuHa=1, nAdj="Y", a=2, 
          alpha1=0.01, beta1=0.15, alpha2=0.1871);
    Data DCSPnormal; Keep Model FSP ESP AveN Power nClassic;
     seedx=1736; seedy=6214; alpha=&alpha; NId=&NId; Nmax=&Nmax; 
     ux=&ux; uy=&uy; sigma=&sigma; model=&Model; N1=&nInterim;
     eSize=abs(&DuHa+NId)/sigma;
     nClassic=round(2*((probit(1-alpha)+probit(1-&beta))/eSize)**2); 
     FSP=0; ESP=0; AveN=0; Power=0;
     Do isim=1 To &nSims;
     nFinal=N1;
       ux1 = Rannor(seedx)*sigma/Sqrt(N1)+ux;
       uy1 = Rannor(seedy)*sigma/Sqrt(N1)+uy;
       T1 = (uy1-ux1+NId)*Sqrt(N1)/2**0.5/sigma;
       p1=1-ProbNorm(T1); 
       If p1>&beta1 then FSP=FSP+1/&nSims;
       If p1<=&alpha1 then do;
          Power=Power+1/&nSims; ESP=ESP+1/&nSims;
       End;
       If p1>&alpha1 and p1<=&beta1 Then Do;
          eRatio = abs(&DuHa/(abs(uy1-ux1)+0.0000001));
          nFinal = min(&Nmax,max(&N0,eRatio**&a*&N0)); 
          If &DuHa*(uy1-ux1+NId) < 0 Then nFinal = N1;
          If &nAdj = "N" then nFinal = &Nmax;
          If nFinal > N1 Then Do;
             ux2 = Rannor(seedx)*sigma/Sqrt(nFinal-N1)+ux ;
             uy2 = Rannor(seedy)*sigma/Sqrt(nFinal-N1)+uy;
             T2 = (uy2-ux2+NId)*Sqrt(nFinal-N1)/2**0.5/sigma;
             p2=1-ProbNorm(T2); 
             If Model="ind" Then TS2=p2; 
             If Model="sum" Then TS2=p1+p2;
             If Model="prd" Then TS2=p1*p2;
             If .<TS2<=&alpha2 Then Power=Power+1/&nSims;
          End;
       End;
       AveN=AveN+nFinal/&nSims;
     End;
     Output;
    run;
    Proc Print Data=DCSPnormal; run; 
    %Mend DCSPnormal;
 
        /*SAS Macro 4.3: Two-Stage Adaptive Design with Survival Endpoint*/
 
    %Macro DCSPSurv(nSims=1000000, model="sum", alpha=0.025, 
          beta=0.2, NId=0, tStd=12, tAcr=4, ux=0, uy=1, DuHa=1,
          nAdj="Y", Nmax=100, N0=100, nInterim=50, a=2, 
          alpha1=0.01, beta1=0.15, alpha2=0.1871);
    Data DCSPSurv; Keep Model FSP ESP AveN Power nClassic;
      seedx=2534; seedy=6762; alpha=&alpha; NId=&NId; 
     Nmax=&Nmax; ux=&ux; uy=&uy; N1=&nInterim; model=&model; 
     Expuxd=exp(-ux*&tStd);  Expuyd=exp(-uy*&tStd);
      sigmax=ux*(1+Expuxd*(1-exp(ux*&tAcr))/(&tAcr*ux))**(-0.5);
      sigmay=uy*(1+Expuyd*(1-exp(uy*&tAcr))/(&tAcr*uy))**(-0.5);
      sigma=((sigmax**2+sigmay**2)/2)**0.5;
      eSize=abs(&DuHa+NId)/sigma;
      nClassic=Round(2*((probit(1-alpha)+Probit(1-&beta))/eSize)**2); 
     FSP=0; ESP=0; AveN=0; Power=0;
      Do isim=1 To &nSims;
      nFinal=N1;
       ux1 = Rannor(seedx)*sigma/Sqrt(N1)+ux;
       uy1 = Rannor(seedy)*sigma/Sqrt(N1)+uy;
       T1 = (uy1-ux1+NId)*Sqrt(N1)/2**0.5/sigma;
       p1=1-ProbNorm(T1); 
       If p1>&beta1 Then FSP=FSP+1/&nSims;
       If p1<=&alpha1 Then do;
          Power=Power+1/&nSims; ESP=ESP+1/&nSims;
       End;
       If p1>&alpha1 and p1<=&beta1 Then Do;
          eRatio=Abs(&DuHa/(Abs(uy1-ux1)+0.0000001));
          nFinal=min(Nmax,max(&N0,eRatio**&a*&N0)); 
          If &DuHa*(uy1-ux1+NId)<0 then nFinal=N1;
          If &nAdj="N" then nFinal=Nmax;
          If nFinal>N1 Then Do;
             ux2 = Rannor(seedx)*sigma/Sqrt(nFinal-N1)+ux ;
             uy2 = Rannor(seedy)*sigma/Sqrt(nFinal-N1)+uy;
             T2 = (uy2-ux2+NId)*Sqrt(nFinal-N1)/2**0.5/sigma;
             p2=1-ProbNorm(T2); 
             If Model="ind" Then TS2=p2; 
             If Model="sum" Then TS2=p1+p2;
             If Model="prd" Then TS2=p1*p2;
             If .<TS2<=&alpha2 Then Power=Power+1/&nSims;
          End;
       End;
       AveN=AveN+nFinal/&nSims;
     End;
     Output;
    Run;
    Proc Print Data=DCSPSurv; Run; 
    %Mend DCSPSurv;
 
        /*SAS Macro 4.4: Event-Based Adaptive Design*/
 
     %Macro DCSPSurv2(nSims=1000000, model="sum", alpha=0.025, 
           beta=0.2, tStd=12, tAcr=4, ux=0.08, uy=0.1, N=100,
           InfoTime=0.5, alpha1=0.01, beta1=0.15, alpha2=0.1871);
     Data DCSPSurv; Keep Model FSP ESP Power AveDs N;
     seed1=2534; seed2=2534; alpha=&alpha; N=&N; ux=&ux; 
       uy=&uy; model=&model; infoTime=&infoTime; tAcr=&tAcr;
     FSP=0; ESP=0; AveDs=0; Power=0; u=(ux+uy)/2;
       Ds=2*N/&tAcr*(&tAcr-(exp(u*&tAcr)-1)/u*exp(-u*&tStd));
       Ds1=Ds*infoTime;
       nFinal=Ds1;
     Do isim=1 To &nSims;
        T1 = Rannor(seed1)+Sqrt(Ds1/2)*log(uy/ux);
        p1=1-ProbNorm(T1); 
        If p1>&beta1 Then FSP=FSP+1/&nSims;
        If p1<=&alpha1 Then do;
           Power=Power+1/&nSims; ESP=ESP+1/&nSims;
        End;
        If p1>&alpha1 and p1<=&beta1 Then Do;
           nFinal=Ds; 
           T2 = Rannor(seed2)+Sqrt((Ds-Ds1)/2)*log(uy/ux);
           p2=1-ProbNorm(T2); 
           If Model="ind" Then TS2=p2; 
           If Model="sum" Then TS2=p1+p2;
           If Model="prd" Then TS2=p1*p2;
           If .<TS2<=&alpha2 Then Power=Power+1/&nSims;
        End;
        AveDs=AveDs+nFinal/&nSims;
     End;
     Output;
     Run;
     Proc Print Data=DCSPSurv; Run; 
     %Mend DCSPSurv2;
 
         /*SAS Macro 4.5: Adaptive Equivalence Trial Design*/
 
    %Macro DCSPEqNormal(nSims=1000000, Model="sum", alpha=0.05, 
    beta=0.2, sigma=0.3, NId=0.2, ux=0, uy=0.1, nInterim=50, Nmax=100, 
    N0=100, DuHa=1, nAdj="Y", a= -2, alpha1=0, beta1=0.2, alpha2=0.3);
     Data DCSPEqNormal; Keep Model FSP ESP AveN Power;
     seedx=1736; seedy=6214; alpha=&alpha; NId=&NId; Nmax=&Nmax; 
     ux=&ux; uy=&uy; sigma=&sigma; model=&Model; N1=&nInterim;
     eSize=abs(&DuHa+NId)/sigma;
     FSP=0; ESP=0; AveN=0; Power=0;
     Do isim=1 To &nSims;
        nFinal=N1;
       ux1 = Rannor(seedx)*sigma/Sqrt(N1)+ux;
        uy1 = Rannor(seedy)*sigma/Sqrt(N1)+uy;
        T11 = (uy1-ux1-NId)*sqrt(N1/2)/sigma; 
        T12 = (uy1-ux1+NId)*sqrt(N1/2)/sigma; 
        p11=Probnorm(T11);      
        p12=Probnorm(-T12);
        p1=max(p11,p12);
            If p1>&beta1 then FSP=FSP+1/&nSims;
        If p1<=&alpha1 then do;
             Power=Power+1/&nSims; ESP=ESP+1/&nSims;
        End;
        If p1>&alpha1 and p1<=&beta1 Then Do;
           eRatio = abs(&DuHa/(abs(uy1-ux1)+0.0000001));
           nFinal = min(&Nmax,max(&N0,eRatio**&a*&N0)); 
           If &DuHa*(uy1-ux1+NId) < 0 Then nFinal = N1;
           If &nAdj = "N" then nFinal = &Nmax;
           If nFinal > N1 Then Do;
              ux2 = Rannor(seedx)*sigma/sqrt(nFinal-N1)+ux ;
              uy2 = Rannor(seedy)*sigma/sqrt(nFinal-N1)+uy;
              T21 = (uy2-ux2-NId)*sqrt(nFinal-N1)/2**0.5/sigma; 
              T22 = (uy2-ux2+NId)*sqrt(nFinal-N1)/2**0.5/sigma;
              p21=Probnorm(T21);
              p22=Probnorm(-T22);      
              p2=max(p21,p22);   
              If Model="ind" Then TS2=p2; 
              If Model="sum" Then TS2=p1+p2;
              If Model="prd" Then TS2=p1*p2;
              If .<TS2<=&alpha2 Then Power=Power+1/&nSims;
           End;
        End;
        AveN=AveN+nFinal/&nSims;
     End;
     Output;
     run;
     Proc Print Data=DCSPEqNormal; run; 
     %Mend DCSPEqNormal;
 
         /*SAS Macro 5.1: Stopping Boundaries with Adaptive Designs*/
    %Macro SB2StgMINP(nSims=10000000, Model="fixedW", w1=0.5, 
          w2=0.5, alpha=0.025, nInterim=50, Nmax=100,
          alpha1=0.01, beta1=0.15, alpha2=0.1871);
    Data SB2StgMINP; Keep Model Power;
    alpha=&alpha; Nmax=&Nmax; Model=&model;
     w1=&w1; w2=&w2; n1=&nInterim; Power=0; seedx=231;
     Do isim=1 To &nSims;
     nFinal=N1;
        T1 = Rannor(seedx);
       p1=1-ProbNorm(T1); 
       If p1<=&alpha1 then do;
          Power=Power+1/&nSims;
       End;
       if p1>&alpha1 and p1<=&beta1 then do;
          T2 = Rannor(seedx);
           If Model^="fixedW" Then do 
             w1=Sqrt(n1/nFinal); 
             w2=Sqrt((1-n1/nFinal)); 
          End;
          Z2=(w1*T1+w2*T2)/Sqrt(w1*w1+w2*w2);
       p2=1-ProbNorm(Z2); 
           If .<p2<=&alpha2 then Power=Power+1/&nSims;
       End;
     End;
    Output;
    Run;
    Proc Print data=SB2StgMINP; Run; 
    %Mend SB2StgMINP;
 
        /*SAS Macro 5.2: Two-Stage Design with Inverse-Normal Method*/
 
    %Macro MINP(nSims=1000000, Model="fixedW", w1=0.5, 
       w2=0.5, alpha=0.025, beta=0.2, sigma=2, NId=0, ux=0, uy=1, 
       nInterim=50, Nmax=100, N0=100, DuHa=1, nAdj="Y", a=2, 
       alpha1=0.01, beta1=0.15, alpha2=0.1871);
    Data MINP; Keep Model FSP ESP AveN Power nClassic PAdj;
      seedx=1736; seedy=6214; alpha=&alpha; NId=&NId; 
     Nmax=&Nmax; Model=&model; w1=&w1; w2=&w2; ux=&ux;
     uy=&uy; sigma=&sigma;  N1=&nInterim; 
     eSize=abs(&DuHa+NId)/sigma;
      nClassic=Round(2*((Probit(1-alpha)+Probit(1-&beta))/eSize)**2); 
      FSP=0; ESP=0; AveN=0; Power=0;
      Do isim=1 To &nSims;
      nFinal=N1;
       ux1 = Rannor(seedx)*sigma/Sqrt(N1)+ux;
       uy1 = Rannor(seedy)*sigma/Sqrt(N1)+uy;
       T1 = (uy1-ux1+NId)*Sqrt(N1)/2**0.5/sigma;
       p1=1-ProbNorm(T1); 
       If p1>&beta1 Then FSP=FSP+1/&nSims;
       If p1<=&alpha1 Then Do;
          Power=Power+1/&nSims; ESP=ESP+1/&nSims;
       End;
       If p1>&alpha1 and p1<=&beta1 Then Do;
          eRatio=Abs(&DuHa/(Abs(uy1-ux1)+0.0000001));
          nFinal=min(&Nmax,max(&N0,eRatio**&a*&N0)); 
          If &DuHa*(uy1-ux1+NId)<0 Then nFinal=N1;
          If &nAdj="N" Then nFinal=&Nmax;
          If nFinal>N1 Then Do;
             ux2 = Rannor(seedx)*sigma/Sqrt(nFinal-N1)+ux ;
             uy2 = Rannor(seedy)*sigma/Sqrt(nFinal-N1)+uy;
             T2 = (uy2-ux2+NId)*Sqrt(nFinal-N1)/2**0.5/sigma;
             If Model^="fixedW" Then Do 
                w1=Sqrt(N1/nFinal); 
                w2=Sqrt((1-N1/nFinal)); 
             End;
             Z2=(w1*T1+w2*T2)/Sqrt(w1*w1+w2*w2);
             p2=1-ProbNorm(Z2); 
             If .<p2<=&alpha2 Then Power=Power+1/&nSims;
          End;
       End;
       AveN=AveN+nFinal/&nSims;
     End;
    PAdj=&alpha1+power-ESP;  ** Stagewise ordering p-value;
     Output;
    Run;
    Proc Print Data=MINP; Run; 
    %Mend MINP;
 
        /*SAS Macro 6.1: N-Stage Adaptive Designs with Normal Endpoint*/
 
    %Macro NStgAdpDsgNor(nSims=1000000, Model="MIP", nStgs=3,
          ux=0, uy=1, NId=0, sigma=2, nAdj="Y", DuHa=1, 
          Nmax=200, N0=150);
    DATA NStgAdpDsg; Set dInput;
    KEEP power Aveux Aveuy AveN FSP1-FSP&nStgs
       ESP1-ESP&nStgs alpha1-alpha&nStgs beta1-beta&nStgs;
    Array Ns{&nStgs}; Array alpha{&nStgs}; Array beta{&nStgs};
    Array ESP{&nStgs}; Array FSP{&nStgs}; 
    seedx=3637; seedy=1624; nStgs=&nStgs; sigma=&sigma; 
    power=0; AveN=0; Aveux=0; Aveuy=0; du=abs(&uy-&ux);
    cumN=0; Do i=1 To nStgs-1; cumN=cumN+Ns{i}; End;
    Do i=1 To nStgs; FSP{i}=0; ESP{i}=0; End;
    Do iSim=1 to &nSims;
     ThisN=0; Thisux=0; Thisuy=0; 
     TS=0; If &Model="MPP" Then TS=1;
     Do i=1 To nStgs;
       uxObs=Rannor(seedx)*sigma/Sqrt(Ns{i})+&ux;    
       uyObs=Rannor(seedy)*sigma/Sqrt(Ns{i})+&uy;    
       Thisux=Thisux+uxObs*Ns{i};
       Thisuy=Thisuy+uyObs*Ns{i};
       ThisN=ThisN+Ns{i};
       Z = (uyObs-uxObs+&NId)*Sqrt(Ns{i}/2)/sigma; 
       pi=1-ProbNorm(Z);
       If &Model="MIP" Then TS=pi;
       If &Model="MSP" Then TS=TS+pi;
       If &Model="MPP" Then TS=TS*pi;
       If TS>beta{i} Then Do; FSP{i}=FSP{i}+1/&nSims;
          Goto Jump; End;
       Else If TS<=alpha{i} then do;
          Power=Power+1/&nSims; ESP{i}=ESP{i}+1/&nSims;
          Goto Jump; End;
       Else If i=1 & &Nadj="Y" Then Do;
         eRatio=&DuHa/(abs(uyObs-uxObs)+0.0000001);
         nFinal=min(&Nmax,max(&N0,eRatio*2*&N0)); 
         If nStgs>1 Then Ns{nStgs}= nFinal-cumN; End;
     End;
    Jump:
    Aveux=Aveux+Thisux/ThisN/&nSims;
    Aveuy=Aveuy+Thisuy/ThisN/&nSims;
    AveN=AveN+ThisN/&nSims;
    End;
    Output;
    Run;
    Proc Print; run;
    %Mend NStgAdpDsgNor;
        
        /*SAS Macro 6.2: N-Stage Adaptive Designs with Binary Endpoint*/
 
    %Macro NStgAdpDsgBin(nSims=1000000, Model="MIP",nStgs=3,
     Px=0, Py=1, NId=0, nAdj="Y", DuHa=1,Nmax=200, N0=150);
    DATA NStgAdpDsg; Set dInput;
    KEEP power AvePx AvePy AveN FSP1-FSP&nStgs 
       ESP1-ESP&nStgs alpha1-alpha&nStgs beta1-beta&nStgs;
    Array Ns{&nStgs}; Array alpha{&nStgs}; Array beta{&nStgs};
    Array ESP{&nStgs}; Array FSP{&nStgs}; 
    seedx=3637; seedy=1624; nStgs=&nStgs; Px=&Px; Py=&Py;
    power=0; AveN=0; AvePx=0; AvePy=0; cumN=0;
    Do i=1 To nStgs-1; cumN=cumN+Ns{i}; End;
    Do i=1 To nStgs; FSP{i}=0; ESP{i}=0; End;
    Do iSim=1 to &nSims;
     ThisN=0; ThisPx=0; ThisPy=0; 
     TS=0; If &Model="MPP" Then TS=1;
    ThisN=0;
     Do i=1 To nStgs;
            PxObs=RanBin(seedx,Ns(i),Px)/Ns(i);
       PyObs=RanBin(seedy,Ns(i),Py)/Ns(i);
       ThisPx=ThisPx+PxObs*Ns(i);
       ThisPy=ThisPy+PyObs*Ns(i);
       ThisN=ThisN+Ns{i};
       sigma=((PxObs*(1-PxObs)+PyObs*(1-PyObs))/2)**0.5;
       Z = (PyObs-PxObs+&NId)*Sqrt(Ns{i}/2)/sigma; 
       pi=1-ProbNorm(Z);
       If &Model="MIP" Then TS=pi;
       If &Model="MSP" Then TS=TS+pi;
       If &Model="MPP" Then TS=TS*pi;
       If TS>beta{i} Then Do; FSP{i}=FSP{i}+1/&nSims; 
                   Goto Jump; End;
       Else If TS<=alpha{i} then do;
          Power=Power+1/&nSims; ESP{i}=ESP{i}+1/&nSims;
                   Goto Jump; End;
       Else If i=1 & &Nadj="Y" Then Do;
          eRatio=&DuHa/(abs(PyObs-PxObs)+0.0000001);
          nFinal=round(min(&Nmax,max(&N0,eRatio*2*&N0))); 
          If nStgs>1 Then Ns{nStgs}= nFinal-cumN; End;
     End;
    Jump:
    AvePx=AvePx+ThisPx/ThisN/&nSims;
    AvePy=AvePy+ThisPy/ThisN/&nSims;
    AveN=AveN+ThisN/&nSims;
    END;
    output;
    RUN;
    proc print; run;
    %Mend NStgAdpDsgBin;
 
        /*SAS Macro 6.3: N-Stage Adaptive Designs with Various Endpoint*/
 
    %Macro TwoArmNStgAdpDsg(nSims=100000, nStgs=3,ux=0, 
       uy=1, NId=0, EP="normal", Model="MSP", Nadj="N",
       DuHa=1, Nmax=300, N0=100, sigma=2, tAcr=10, tStd=24);
    DATA NStgAdpDsg; Set dInput;
    KEEP power Aveux Aveuy AveN FSP1-FSP&nStgs
       ESP1-ESP&nStgs alpha1-alpha&nStgs beta1-beta&nStgs;
    Array Ns{&nStgs}; Array alpha{&nStgs}; Array beta{&nStgs};
    Array ESP{&nStgs}; Array FSP{&nStgs}; Array Ws{&nStgs};
    Array sumWs{&nStgs}; Array TSc{&nStgs};
    seedx=3637; seedy=1624; Model=&Model; nStgs=&nStgs;
    sigma=&sigma; power=0; AveN=0; Aveux=0; Aveuy=0;
    cumN=0; Do i=1 To nStgs-1; cumN=cumN+Ns{i}; End;
    Do k=1 To nStgs;
        sumWs{k}=0; Do i=1 To k; 
                   sumWs{k}=sumWs{k}+Ws{i}**2; End; 
       sumWs{k}=Sqrt(sumWs{k});
    End;
    * Calcate the standard deviation, sigma for different endpoints *;
     u=(&ux+&uy)/2;
     if &EP="normal" Then sigma=&sigma;
     if &EP="binary" Then sigma=(u*(1-u))**0.5;
     if &EP="survival" Then Do;
       expud=exp(-u*&tStd);
       sigma=u*(1+expud*(1-exp(u*&tAcr))/(&tAcr*u))**(-0.5);
     End;
    Do i=1 To nStgs; FSP{i}=0; ESP{i}=0; End;
    Do iSim=1 to &nSims;
     ThisN=0; Thisux=0; Thisuy=0; 
     Do i=1 To nStgs; TSc{i}=0; End; 
     TS=0; If &Model="MPP" Then TS=1;
     Do i=1 To nStgs;
       uxObs=Rannor(seedx)*sigma/Sqrt(Ns{i})+&ux;    
       uyObs=Rannor(seedy)*sigma/Sqrt(Ns{i})+&uy;    
       Thisux=Thisux+uxObs*Ns{i};
       Thisuy=Thisuy+uyObs*Ns{i};
       ThisN=ThisN+Ns{i};
       TS0 = (uyObs-uxObs+&NId)*Sqrt(Ns{i}/2)/sigma; 
    If Model="MIP" Then TS=1-ProbNorm(TS0);
    If Model="MSP" Then TS=TS+(1-ProbNorm(TS0));
    If Model="MPP" Then TS=TS*(1-ProbNorm(TS0));
    If Model="WZ" Then Do;
       Do k=i to nStgs; 
          TSc{k}=TSc{k}+Ws{i}/sumWs{k}*TS0;
       End;
       TS=1-ProbNorm(TSc{i}); 
    End;
    If Model="UWZ" Then Do; 
       TS0=((Thisuy-Thisux)/ThisN+&NId)*Sqrt(ThisN/2)/sigma;
       TS=1-ProbNorm(TS0); 
    End;
       If TS>beta{i} Then Do; FSP{i}=FSP{i}+1/&nSims; 
                      Goto Jump; End;
       Else If TS<=alpha{i} then do;
          Power=Power+1/&nSims; ESP{i}=ESP{i}+1/&nSims;
                      Goto Jump; End;
       Else If i=1 & &Nadj="Y" Then Do;
          eRatio=&DuHa/(abs(uyObs-uxObs)+0.0000001);
          nFinal=min(&Nmax,max(&N0,eRatio*2*&N0)); 
          If nStgs>1 Then Ns{nStgs}= nFinal-cumN; End;
     End;
    Jump:
    Aveux=Aveux+Thisux/ThisN/&nSims;
    Aveuy=Aveuy+Thisuy/ThisN/&nSims;
    AveN=AveN+ThisN/&nSims;
    End;
    Output;
    Run;
    Proc Print; Run;
    %Mend TwoArmNStgAdpDsg;
 
        /*SAS Macro 7.1: Conditional Power*/
 
    %Macro ConPower(EP="normal", Model="MSP", alpha2=0.205, 
     ux=0.2, uy=0.4, sigma=1, n2=100, p1=0.8, w1=1, w2=1);
    ** cPower=Two stage conditional power. eSize=delta/sigma;
    data cPower;
     a2=&alpha2; Model=&Model;
     u=(&ux+&uy)/2;
     w1=&w1/sqrt(&w1**2+&w2**2);
     w2=&w2/sqrt(&w1**2+&w2**2);
     If &EP="normal" Then sigma=&sigma;
     If &EP="binary" Then sigma=(u*(1-u))**0.5;
     eSize=(&uy-&ux)/sigma;
    If Model="MIP" Then BFun=Probit(1-a2);
    If Model="MSP" Then BFun=Probit(1-max(0.0000001,a2-&p1)); 
    If Model="MPP" Then BFun=Probit(1-a2/&p1);
    If Model="LW" Then BFun=(Probit(1-a2)- w1*Probit(1-&p1))/w2;
    cPower=1-ProbNorm(BFun-eSize*sqrt(&n2/2)); 
    Run;
    Proc Print data=cPower; Run;
    %Mend ConPower;
 
        /*SAS Macro 7.2: Sample-Size Based on Conditional Power*/
 
    %Macro nByCPower(nAdjModel, alpha2, eSize, 
        cPower, p1, w1, w2, n2New);
     a2=&alpha2; 
     If &nAdjModel="MIP" Then BFun=Probit(1-a2);
     If &nAdjModel="MSP" Then BFun=Probit(1-max(0.0000001,a2-&p1)); 
     If &nAdjModel="MPP" Then BFun=Probit(1-a2/&p1);
     If &nAdjModel="LW" Then 
       BFun=(Probit(1-a2)- &w1*Probit(1-&p1))/&w2;
     &n2New=2*((BFun-Probit(1-&cPower))/&eSize)**2; *n per group;
     %Mend nByCPower;
 
         /*SAS Macro 9.1: Adaptive Design with Sample-Size Re-estimation */
 
    %Macro nByCPowerUB(Model, a2, eSize, cPower, p1, w1, w2, n2New);
     If &Model="MIP" Then BFun=Probit(1-&a2);
      If &Model="MSP" Then BFun=Probit(1-max(0.0000001,&a2-&p1)); 
     If &Model="MPP" Then BFun=Probit(1-&a2/&p1);
      If &Model="LW" Then 
       BFun=(Probit(1-&a2)- &w1*Probit(1-&p1))/&w2;
     &n2New=2*((BFun-Probit(1-&cPower))/&eSize)**2; *n for group x;
    %Mend nByCPowerUB;
    
    %Macro TwoArmNStgAdpDsg(nSims=1000000, nStgs=2,ux=0, uy=1, 
     NId=0, NItype="FIXED", EP="normal", Model="MSP", 
     Nadj="Y", cPower=0.9, DuHa=1, Nmax=300, N2min = 150, 
     nMinIcr=1, sigma=3, tAcr=10, tStd=24, nRatio=1);
    DATA NStgAdpDsg; Set dInput;
    KEEP Model power Aveux Aveuy AveTotalN FSP1-FSP&nStgs 
       ESP1-ESP&nStgs;
    Array Ns{&nStgs}; Array alpha{&nStgs}; Array beta{&nStgs};
    Array ESP{&nStgs}; Array FSP{&nStgs}; Array Ws{&nStgs};
    Array sumWs{&nStgs}; Array TSc{&nStgs};
    Model=&Model; cPower=&cPower; nRatio=&nRatio; NId=&NId; 
    NItype=&NItype; N2min=&N2min; nStgs=&nStgs; sigma=&sigma; 
    power=0; AveTotalN=0; Aveux=0; Aveuy=0; ux=&ux; uy=&uy;
    cumN=0; Do i=1 To nStgs-1; cumN=cumN+Ns{i}; End;
    N0=CumN+Ns{nStgs};
    Do k=1 To nStgs;
        sumWs{k}=0; Do i=1 To k; sumWs{k}=sumWs{k}+Ws{i}**2; End; 
        sumWs{k}=sqrt(sumWs{k});
    End;
    * Calcate the standard deviation, sigma for different endpoints *;
     If &EP="normal" Then Do sigmax=&sigma; sigmay=&sigma; End;
     If &EP="binary" Then Do 
          sigmax=Sqrt(&ux*(1-&ux)); sigmay=Sqrt(&uy*(1-&uy)); 
     End;
     If &EP="survival" Then Do
         sigmax=ux*sqrt(1+exp(-ux*&tStd)*(1-exp(ux*&tAcr))/(&tAcr*ux));
         sigmay=uy*sqrt(1+exp(-uy*&tStd)*(1-exp(uy*&tAcr))/(&tAcr*uy));
     End;
     If NItype="PCT" Then sigmax=(1-NId)*sigmax;
         Do i=1 To nStgs; FSP{i}=0; ESP{i}=0; End;
    Do iSim=1 to &nSims;
     ThisNx=0; ThisNy=0; Thisux=0; Thisuy=0;
     Do i=1 To nStgs; TSc{i}=0; End; 
     TS=0; If &Model="MPP" Then TS=1;
          Do i=1 To nStgs;
        uxObs=Rannor(746)*sigmax/sqrt(Ns{i})+&ux;    
        uyObs=Rannor(874)*sigmay/sqrt(nRatio*Ns{i})+&uy;    
                  Thisux=Thisux+uxObs*Ns{i};
        Thisuy=Thisuy+uyObs*nRatio*Ns{i};
        If NItype="PCT" Then NId=uxObs*&NId;
        ThisNx=ThisNx+Ns{i};
        ThisNy=ThisNy+Ns{i}*nRatio;
        StdErr=(sigmax**2/Ns{i}+sigmay**2/(nRatio*Ns{i}))**0.5;   
        TS0 = (uyObs-uxObs+NId)/StdErr; 
                  If Model="MIP" Then TS=1-ProbNorm(TS0);
        If Model="MSP" Then TS=TS+(1-ProbNorm(TS0));
        If Model="MPP" Then TS=TS*(1-ProbNorm(TS0));
        If Model="LW" Then Do; 
       Do k=i to nStgs; TSc{k}=TSc{k}+Ws{i}/sumWs{k}*TS0; End;
       TS=1-ProbNorm(TSc{i}); 
        End;
        If Model="UWZ" Then Do; 
       StdErr=(sigmax**2/ThisNx+sigmay**2/ThisNy)**0.5;   
       TS0=(Thisuy/ThisNy-Thisux/ThisNx+NId)/StdErr;
       TS=1-ProbNorm(TS0); 
        End;
        If TS>beta{i} Then Do; FSP{i}=FSP{i}+1/&nSims; Goto Jump; End;
        Else If TS<=alpha{i} then do;
       Power=Power+1/&nSims; ESP{i}=ESP{i}+1/&nSims; 
          Goto Jump; End;
        Else If nStgs>1 & i=1 & &Nadj="Y" Then Do;
       eSize=&DuHa/(abs(uyObs-uxObs)+0.0000001);
        nFinal=min(&Nmax, max(N0,eSize*Abs(eSize)*N0));
                 If nStgs=2 Then do;
       * eSize=(uyObs-uxObs+NId)/((sigmax+sigmay)*sqrt(0.5+0.5/nRatio)); 
        eSize=(uyObs-uxObs+NId)/Sqrt(0.5*sigmax**2+0.5*sigmay**2/nRatio);
       %nByCPowerUB(Model, alpha{2}, eSize, cPower, TS, 
          ws{1}, ws{2}, n2New);
       nFinal=Round(min(&Nmax,ns{1}+n2New+&nMinIcr/2), &nMinIcr); 
       nFinal=max(N2min+Ns{1},nFinal); 
       End;
              If nStgs>1 Then Ns{nStgs}= max(1,nFinal-cumN); 
        End;
     End;
         Jump:
    Aveux=Aveux+Thisux/ThisNx/&nSims;
    Aveuy=Aveuy+Thisuy/ThisNy/&nSims;
    AveTotalN=AveTotalN+(ThisNx+ThisNy)/&nSims; 
    End;
    Output;
    Run;
    Proc Print; Run;
    %Mend TwoArmNStgAdpDsg;
 
        /*SAS Macro 11.1: Two-Stage Drop-Loser Adaptive Design*/
 
    %Macro DrpLsrNRst(nSims=100000, CntlType="strong", nArms=5, 
       alpha=0.025, beta=0.2, NId=0, cPower=0.9, nInterim=50,
       Nmax=150, nAdj="Y", alpha1=0.01, beta1=0.15, 
       alpha2=0.1871, EP="normal", sigma=1, tStd=24, tAcr=10);
    Data DrpLsrNRst; Set dInput; 
     Keep FSP ESP AveN Power cPower Nmax;
     Array us{&nArms}; Array u1{&nArms}; 
     Array u2{&nArms}; Array cs{&nArms};
     seedx=1736; seedy=6214; alpha=&alpha; NId=&NId;
     Nmax=&Nmax; nArms=&nArms; n1=&nInterim; sigma=&sigma; 
     cPower=&cPower; FSP=0; ESP=0; AveN=0; Power=0;
     If &EP="mean" Then sigma=&sigma;
     If &EP="binary" Then sigma=(us1*(1-us{1}))**0.5;
     If &EP="survival" Then 
       sigma=us{1}*(1+exp(-us{1}*&tStd)*(1-exp(us{1}*&tAcr))
                /(&tAcr*us{1}))**(-0.5);
          Do isim=1 to &nSims;
     TotalN=nArms*n1; uMax=us{1}; Cntrst=0; SumSqc=0;
       Do i=1 To nArms;
       u1{i}=Rannor(seedx)*sigma/Sqrt(n1)+us{i};
       If u1{i}>uMax Then Do uMax=u1{i}; iMax=i; End;
       Cntrst=Cntrst+cs{i}*u1{i};
       SumSqc=SumSqc+cs{i}*cs{i};
        End;
       Z1 = Cntrst*Sqrt(n1)/Sqrt(SumSqc)/sigma;
       p1=1-ProbNorm(Z1); 
       * Bonferroni Adjustment;
       If &CntlType="strong" Then 
         p1=(nArms-1)*(1-ProbNorm((uMax-us{1})/sigma*Sqrt(n1/2)));
       If p1>&beta1 Then FSP=FSP+1/&nSims; 
       If p1<=&alpha1 Then do;
       Power=Power+1/&nSims; ESP=ESP+1/&nSims;
       End;
       If iMax=1 Then Goto myJump;
       If p1>&alpha1 and p1<=&beta1 Then do;
       BF=Probit(1-max(0,&alpha2-p1))-Probit(1-cPower);
       n2=2*(sigma/(u1{iMax}-u1{1})*BF)**2;
       nFinal=min(n1+n2, Nmax);
       If &nAdj="N" Then nFinal=Nmax;
       If nFinal>n1 Then do;
          TotalN=2*(nFinal-n1)+nArms*n1;
          u2{1}=Rannor(seedx)*sigma/Sqrt(nFinal-n1)+us{1};
           u2{iMax}=Rannor(seedy)*sigma/Sqrt(nFinal-n1)+us{iMax};
           T2=(u2{iMax}-u2{1}+NId)*Sqrt(nFinal-n1)/2**0.5/sigma;
       p2=1-ProbNorm(T2); TS2=p1+p2;
           If .<TS2<=&alpha2 Then Power=Power+1/&nSims;
       End;
       End;
    myJump:
       AveN=AveN+TotalN/&nSims;
     End;
     Output;
    Run;
    Proc Print Data=DrpLsrNRst (obs=1); Run; 
    %Mend DrpLsrNRst;
 
        /*SAS Macro 12.1: Biomarker-Adaptive Design*/
 
    %Macro BMAD(nSims=100000, cntlType="strong", nStages=2, 
       u0p=0.2, u0n=0.1, sigma=1, np1=50, np2=50, nn1=100,
       nn2=100, alpha1=0.01, beta1=0.15,alpha2=0.1871);
    Data BMAD; 
     Keep FSP ESP Power AveN pPower oPower;
     seedx=1736; seedy=6214; u0p=&u0p; u0n=&u0n; np1=&np1;
       np2=&np2; nn1=&nn1; nn2=&nn2; sigma=&sigma; 
       FSP=0; ESP=0;Power=0; AveN=0; pPower=0; oPower=0;
     Do isim=1 to &nSims;
        up1=Rannor(seedx)*sigma/Sqrt(np1)+u0p;
        un1=Rannor(seedy)*sigma/Sqrt(nn1)+u0n;
        uo1=(up1*np1+un1*nn1)/(np1+nn1);
        Tp1=up1*np1**0.5/sigma; To1=uo1*(np1+nn1)**0.5/sigma;
        T1=Max(Tp1,To1); p1=1-ProbNorm(T1); 
        If &cntlType="strong" Then p1=2*p1; *Bonferroni;
        If p1>&beta1 Then FSP=FSP+1/&nSims; 
        If p1<=&alpha1 Then Do;
            Power=Power+1/&nSims; ESP=ESP+1/&nSims;
            If Tp1>To1 Then pPower=pPower+1/&nSims;
            If Tp1<=To1 Then oPower=oPower+1/&nSims;
       End;
       AveN=AveN+2*(np1+nn1)/&nSims;
       If &nStages=2 And p1>&alpha1 And p1<=&beta1 Then Do;
           up2=Rannor(seedx)*sigma/Sqrt(np2)+u0p;
           un2=Rannor(seedy)*sigma/Sqrt(nn2)+u0n;
           uo2=(up2*np2+un2*nn2)/(np2+nn2);
           Tp2=up2*np2**0.5/sigma;  To2=uo2*(np2+nn2)**0.5/sigma;
           If Tp1>To1 Then Do; 
          T2=Tp2; AveN=AveN+2*np2/&nSims;
         End;
           If Tp1<=To1 Then Do;
          T2=To2; AveN=AveN+2*(np2+nn2)/&nSims;
        End;
           p2=1-ProbNorm(T2); Ts=p1+p2;
           If .<TS<=&alpha2 Then Do;
          Power=Power+1/&nSims;
          If Tp1>To1 Then pPower=pPower+1/&nSims;
          If Tp1<=To1 Then oPower=oPower+1/&nSims;
          End;
       End;
     End;
    Run;
    Proc Print Data=BMAD (obs=1); Run; 
    %Mend BMAD;
 
        /*SAS Macro 14.1: Randomized Play-the-Winner Design*/
 
    %Macro RPW(nSims=100000, Zc=1.96, nSbjs=200, nAnlys=3,
       RR1=0.2, RR2=0.3, a0=1, b0=1, a1=1, b1=1, nMin=1);
    Data RPW; Keep nSbjs aveP1 aveP2 Zc Power;
    seed1=364; seed2=894; Power=0; aveP1=0; aveP2=0;   
    Do isim=1 to &nSims;
     nResp1=0; nResp2=0; a0=&a0; b0=&b0; Zc=&Zc; N1=0; N2=0;
     nSbjs=&nSbjs; nAnlys=&nAnlys; nMax=nSbjs-&nMin;
     a=a0; b=b0; r0=a/(a+b);
     Do iSbj=1 To nSbjs;
        If Mod(iSbj,Round(nSbjs/nAnlys))=0 Then r0=a/(a+b);
       rnAss=Ranuni(seed1);
       If (rnAss < r0 And N1<nMax) Or N2>=nMax Then 
          Do;
             N1=N1+1; rnRep=Ranuni(seed2);
             if rnRep <=&RR1 Then Do;
                nResp1=nResp1+1; a=a+&a1;
             End;
          End;
       Else 
          Do;
             N2=N2+1; rnRep=Ranuni(seed2);
             If rnRep <=&RR2 Then Do;
                nResp2=nResp2+1; b=b+&b1;
             End;
          End;
     End;
     p1=nResp1/N1; p2=nResp2/N2;
     aveP1=aveP1+p1/&nSims; aveP2=aveP2+p2/&nSims;
     sigma1=sqrt(p1*(1-p1)); sigma2=sqrt(p2*(1-p2));
     Sumscf=sigma1**2/(N1/(N1+N2))+sigma2**2/(N2/(N1+N2));
     TS = (p2-p1)*Sqrt((N1+N2)/sumscf); 
     If TS>Zc Then Power=Power+1/&nSims; 
    End;
    Output;
    Run;
    Proc Print data=RPW; Run;
    %Mend RPW;
 
        /*SAS Macro 14.2: Binary Response-Adaptive Randomization*/
 
    %Macro RARBin(nSims=1000, nPts=200, nArms=5, b=1, m=1, Zc=1.96);
     Data RARBin; Set DataIn;
     Keep nPts AveN1-AveN&nArms AveU1-AveU&nArms PowerMax;
     Array Ns{&nArms}; Array uObs{&nArms}; Array rP{&nArms}; 
     Array nRsps{&nArms}; Array a0{&nArms}; Array CrP{&nArms}; 
     Array us{&nArms}; Array AveU{&nArms}; Array AveN{&nArms}; 
     PowerMax=0; nArms=&nArms; nPts=&nPts; 
     Do i=1 To nArms; AveU{i}=0; AveN{i}=0; End; 
     Do isim=1 to &nSims;
         Do i=1 To nArms; nRsps{i}=0; uObs{i}=0; Ns{i}=0; Crp{i}=0; End;
         Do iSubject=1 to nPts;
        Do i=1 To nArms; rP{i}=a0{i}+&b*uObs{i}**&m; End;
        Suma=0; Do i=1 To nArms; Suma=Suma+rP{i}; End;
        Do i=1 To nArms; rP{i}=rP{i}/Suma; End;
        Do iArm=1 To nArms; CrP{iArm}=0;
           Do i=1 To iArm; CrP{iArm}=CrP{iArm}+rP{i}; End;
        End;
        rn=ranuni(5236); cArm=1;
        Do iArm=2 To nArms; 
           IF CrP{iArm-1}<rn<CrP{iArm} Then cArm=iArm; 
        End; 
        Ns(cArm)= Ns(cArm)+1; 
        * For Binary response;
        If ranuni(8364)<us{cArm} Then nRsps{cArm}=nRsps{cArm}+1;
        Do i=1 To nArms; uObs{i}=nRsps{i}/max(Ns{i},1); End;
         End;
       uMax=uObs{1}; 
         Do i=1 to &nArms; If uObs{i}>=uMax Then iMax=i; End;
         Se2=0;
         Do i=1 to &nArms; 
       Se2=Se2+uObs{i}*(1-uObs{i})/max(Ns{i},1)*2/nArms; 
         End;
         TSmax=(uObs{iMax}-uObs{1})*(Ns(1)+Ns{iMax})/2
                      /(nPts/nArms)/Se2**0.5; 
         If TSmax>&Zc then PowerMax=PowerMax+1/&nSims; 
         Do i=1 To nArms;
        AveU{i}=AveU{i}+uObs{i}/&nSims; 
        AveN{i}=AveN{i}+Ns{i}/&nSims; 
        End;
     End;
     Output;
     Run;
     Proc Print Data=RARBin; Run; 
     %Mend RARBin;
 
         /*SAS Macro 14.3: Normal Response-Adaptive Randomization*/
 
    %Macro RARNor(nSims=100000, nPts=100, nArms=5, b=1, m=1,
                             CrtMax=1.96);
    Data RARNor; Set DataIn; 
    Keep nPts AveN1-AveN&nArms AveU1-AveU&nArms PowerMax;
    Array Ns{&nArms}; Array AveN{&nArms}; Array uObs{&nArms}; 
    Array rP{&nArms}; Array AveU{&nArms}; Array cuObs{&nArms}; 
    Array a0{&nArms}; Array CrP{&nArms}; 
    Array us{&nArms}; Array s{&nArms}; 
    PowerMax=0; nArms=&nArms; nPts=&nPts;
    Do i=1 To nArms; AveU{i}=0; AveN{i}=0; End;   
    Do isim=1 to &nSims; 
     Do i=1 To nArms; cuObs{i}=0; uObs{i}=0; Ns{i}=0; Crp{i}=0; End;
     Do iSubject=1 to nPts;
        Do i=1 To nArms; rP{i}=a0{i}+&b*uObs{i}**&m; End;
        Suma=0; Do i=1 To nArms; Suma=Suma+rP{i}; End;
        Do i=1 To nArms; rP{i}=rP{i}/Suma; End;
       Do iArm=1 To nArms; CrP{iArm}=0;
          Do i=1 To iArm; CrP{iArm}=CrP{iArm}+rP{i}; End;
       End;
       rn=ranuni(5361); cArm=1;
     Do iArm=2 To nArms; 
          IF CrP{iArm-1}<rn<CrP{iArm} Then cArm=iArm; 
     End; 
       Ns(cArm)= Ns(cArm)+1; 
     * For Normal response;
       u=Rannor(361)*s{cArm}+us{cArm};
     cuObs{cArm}=cuObs{cArm}+u; 
       Do i=1 To nArms; uObs{i}=cuObs{i}/max(Ns{i},1); End;
     End;
     se2=0;
     * Assume sigma unknown for simplicity;
     Do i=1 To nArms; se2=se2+s{i}**2/max(Ns{i},1)*2/nArms; End; 
     uMax=uObs{1}; 
     Do i=1 To &nArms; If uObs{i}>=uMax Then iMax=i; End;
     TSmax=(uObs{iMax}-uObs{1})*(Ns(1)+Ns{iMax})/2
                   /(nPts/nArms)/se2**0.5; 
     If TSmax>&CrtMax then PowerMax=PowerMax+1/&nSims; 
     Do i=1 To nArms;
        AveU{i}=AveU{i}+uObs{i}/&nSims; 
        AveN{i}=AveN{i}+Ns{i}/&nSims; 
     End;
    End;
    Output;
    Run;
    Proc Print data=RARNor; Run; 
    %Mend RARNor;
 
        /*SAS Macro 15.1: 3+3 Dose-Escalation Design*/
 

   %Macro TER3p3(nSims=50000, nLevels=10);

Data TER; Set dInput; Keep AveMTD AveNPts AveNRsps;

Array nPts{&nLevels}; Array nRsps{&nLevels}; Array RspRates{&nLevels};

AveMTD=0; AveNPts=0; AveNRsps=0;  nLevels=&nLevels;

Do iSim=1 to &nSims;

  Do i=1 To nLevels; nPts{i}=0; nRsps{i}=0; End;

  seedn=Round((Ranuni(281)*100000000));

  iLevel=1; TotPts=0; TotRsps=0;

Looper:

  If iLevel>&nLevels | nPts{iLevel}=6 Then Goto Finisher; 

   nPts{iLevel}=nPts{iLevel}+3;

  rspRate=RspRates{iLevel};

  Rsp=RANBIN(seedn,3,rspRate);

  nRsps{iLevel}=nRsps{iLevel}+Rsp;

  TotPts=TotPts+3; TotRsps=TotRsps+Rsp;

  If nPts(iLevel)=3 & nRsps{iLevel}=0 Then Do;

      iLevel=iLevel+1; Goto Looper;

  End;

  If nPts(iLevel)=3 & nRsps{iLevel}=1 Then Goto Looper;

  If nPts(iLevel)=6 & nRsps{iLevel}<=1 Then Do;

      iLevel=iLevel+1; Goto Looper;

  End;

Finisher:

  MTD=Min(iLevel-1, nLevels);

  AveMTD=AveMTD+MTD/&nSims;

  AveNPts=AveNPts+totPts/&nSims;

  AveNRsps=AveNRsps+TotRsps/&nSims;

End;

Output;

Run;

Proc Print Data=TER; Run;

%Mend TER3p3;

 
 
        /*SAS Macro 15.2: Continual Reassessment Method*/
 
    %Macro CRM(nSims=100, nPts=30, nLevels=10, b=100, 
       aMin=0.1, aMax=0.3, MTRate=0.3, nIntPts=100);
    Data CRM; Set DInput; Keep nPts nLevels AveMTD SdMTD DLTS; 
    Array nPtsAt{&nLevels}; Array nRsps{&nLevels}; Array g{&nIntPts};
    Array Doses{&nLevels}; Array RRo{&nLevels}; Array RR{&nLevels}; 
    Array g0{&nIntPts};  
    seed=2736; nLevels=&nLevels; nPts=&nPts; DLTs=0; 
    AveMTD=0;VarMTD=0; dx=(&aMax-&aMin)/&nIntPts;
    Do k=1 To &nIntPts; g0{k}=g{k}; End; 
    Do iSim=1 to &nSims;
    Do k=1 To &nIntPts; g{k}=g0{k}; End; 
     Do i=1 To nLevels; nPtsAt{i}=0; nRsps{i}=0; End;
     iLevel=1; 
     Do iPtient=1 To nPts;
     iLevel=min(iLevel, &nLevels); Rate=RRo{iLevel};
     nPtsAt{iLevel}=nPtsAt{iLevel}+1;
       r=Ranbin(seed,1,Rate); nRsps{iLevel}=nRsps{iLevel}+r;
    ** Posterior distribution of a;
       c=0;
       Do k=1 To &nIntPts;
          ak=&aMin+k*dx; 
          Rate=1/(1+&b*Exp(-ak*doses{iLevel}));
          If r>0 Then L=Rate; Else L=(1-Rate);
          g{k}=L*g{k}; c=c+g{k}*dx;
       End;
       Do k=1 to &nIntPts; g{k}=g{k}/c; End;
    ** Predict response rate and current MTD;
        MTD=iLevel; MinDR=1;
        Do i=1 To nLevels;
           RR{i}=0;
          Do k=1 To &nIntPts;
          ak=&aMin+k*dx;
        RR{i}= RR{i}+1/(1+&b*Exp(-ak*doses{i}))*g{k}*dx;
          End;
          DR=Abs(&MTRate-RR{i});
          If .<DR <MinDR Then 
             Do; MinDR = DR; iLevel=i; MTD=i; End;
       End;
     End;
     Do i=1 To nLevels; 
       DLTs=DLTs+nRsps{i}/&nSims; 
     End;
     AveMTD=AveMTD+MTD/&nSims;
     VarMTD=VarMTD+MTD**2/&nSims;
    End;
    SdMTD=(VarMTD-AveMTD**2)**0.5;
    Output;
    Run;
    Proc Print Data=CRM; run;
    %Mend CRM;
 
        /*SAS Macro 16.1: Simon Two-Stage Futility Design*/
 
    %Macro TwoStageDesign(n=50, po=0.15, pa=0.3, n1=20, n1c=2, 
    alpha0=0.1);
    Data TwoStageBin;
     retain alpha power;
     drop i p1o p2o p1a p2a n2;
     n1=&n1; n1c=&n1c; po=&po; pa=&pa; * Remove "&".;
     do n=round(0.7*&n) to round(1.5*&n);
        n2=n-&n1; 
       do nc=n1c to n;
          alpha=0;
          power=0;
          do i=max(n1c,nc-n2) to n1;
              p1o=ProbBnml(po, n1,i)-ProbBnml(po, n1,i-1);
             p2o=1;
             if nc-i>0 then p2o=1-ProbBnml(po, n2,nc-i-1);
             alpha=alpha+p1o*p2o;
             p1a=ProbBnml(pa, n1,i)-ProbBnml(pa, n1,i-1);
             p2a=1;
             if nc-i>0 then p2a=1-ProbBnml(pa, n2,nc-i-1);
             power=power+p1a*p2a;
          end;
          if alpha>0.8*&alpha0 && alpha<1.2*&alpha0 then do;
             PrEFSHo=ProbBnml(po, n1,n1c-1);
             ExpNo=PrEFSHo*n1+(1-PrEFSHo)*n;
             PrEFSHa=ProbBnml(pa, n1,n1c-1);
             ExpNa=PrEFSHa*n1+(1-PrEFSHa)*n;
             output;
          end;
       end;
     end;
    run;
    proc print; run;
    run;
    %Mend TwoStageDesign;
 
 

 

/* SAS Macro A.1 Mixed Exponential Distribution */

%Macro RanVars(nObs=100, Lamda1=1, Lamda2=1.5, w1=0.6, w2=0.4);

Data RVars;

Drop iObs;

Do iObs=1 To &nObs;

xMixEXP=&w1/&Lamda1*Ranexp(782)+&w2/&Lamda2*Ranexp(323);

Output;

End;

Run;

Proc Print Data=RVars; Run;

%Mend RanVars;

 
 

            /* SAS Macro A.2 Multi-Variate Normal Distribution */

%Macro RanVarMNor(sSize=4, nVars=2, nObs=10);

Data nVars; SET CorrMatrix; keep x1-x&nVars;

Array a{&sSize}; Array xNor{&nVars}; Array x{&nVars};

Array ss{&sSize}; * Correlation matrix;

* Checovsky decomposition;

Do k=1 To &nVars; Do L=1 To k;

Saa=0; Sa2=0;

Do j=1 to L-1;

Saa=Saa+a{&nVars*(k-1)+j}*a{&nVars*(L-1)+j};

Sa2=Sa2+a{&nVars*(L-1)+j}*a{&nVars*(L-1)+j};

End;

nkL=&nVars*(k-1)+L;

a{nkL}=(ss{nkL}-Saa)/Sqrt(ss{&nVars*(L-1)+L}-Sa2);

End; End;

Do iObs=1 to &nObs;

Do iVar=1 to &nVars; xNOR{iVar}=Rannor(762); End;

Do iVar=1 to &nVars;

x{iVar}=0;

Do i=1 To iVar;

x{iVar}=x{iVar}+a{&nVars*(iVar-1)+i}*xNor{i};

End;

End;

Output;

End;

Run;

Proc Print data=nVars(obs=100); Run;

Proc corr data=nVars; Run;

%Mend RanVarMNor;