data step1;

** You may change the number of simulations from 1000000 to 10000;

** if you wish to reduce the cpu time;
nsimul=1000000;

** You may change alpha to the level say 0.05 as you wish;

** You may change za to the percentile say 1.645 as you wish;
alpha=0.025; za=1.96; zb=probit(0.8);

** If you wish to change the range for the following pxyz;

** Make sure that all pxyz fall between 0 and 1;

** Otherwise there will an error;
do lowbdor=0.50; ** lower bound for non-inferiority of the odds ratio;
do or=0.90,1.0,1.1;  ** the underlying OR;
do p011=0.05,0.10,0.20;
do p101=0.15,0.25;
   pc1=1-p011-p101;  ** probability of concordance in group 1;
do p012=0.10;
   p102=or**2*p101*p012/p011;
   pc2=1-p012-p102;
   vlor=(1/p011+1/p101+1/p012+1/p102)/4;
  aprxn=int((za+zb)**2*vlor/(log(lowbdor)-log(or))**2)+1;
   n=aprxn;
   sumrej=0;
   sumrej1=0;
   p102=lowbdor**2*p101*p012/p011;
   pc2=1-p012-p102;
    do i=1 to nsimul;
     sn011=0; sn101=0; sn012=0; sn102=0;
     do nx=1 to n;
     x=rantbl(1432151,p011,p101,pc1);
     y=rantbl(1431171,p012,p102,pc2);
     if x eq 1 then sn011+1;
     if x eq 2 then sn101+1;
     if y eq 1 then sn012+1;
     if y eq 2 then sn102+1;
     end;  ** end of do nx=1 to n;
     ep011=(sn011+0.5)/n;
     ep101=(sn101+0.5)/n;
     ep012=(sn012+0.5)/n;
     ep102=(sn102+0.5)/n;
     estor=sqrt((ep011*ep102)/(ep101*ep012));
    ztest=(log(estor)-log(lowbdor))/(sqrt(1/(sn011+0.5)+1/(sn101+0.5)+
                                          1/(sn012+0.5)+1/(sn102+0.5))/2);
    if ztest ge 1.96 then sumrej1+1;
    end;  ** end of do nsimul;
   err1=sumrej1/nsimul;
   aprxerr1=err1;
   p102=or**2*p101*p012/p011;
   pc2=1-p012-p102;
    do i=1 to nsimul;
     sn011=0; sn101=0; sn012=0; sn102=0;
     do nx=1 to n;
     x=rantbl(1432151,p011,p101,pc1);
     y=rantbl(1431171,p012,p102,pc2);
     if x eq 1 then sn011+1;
     if x eq 2 then sn101+1;
     if y eq 1 then sn012+1;
     if y eq 2 then sn102+1;
     end;  ** end of do nx=1 to n;
     ep011=(sn011+0.5)/n;
     ep101=(sn101+0.5)/n;
     ep012=(sn012+0.5)/n;
     ep102=(sn102+0.5)/n;
     estor=sqrt((ep011*ep102)/(ep101*ep012));
    ztest=(log(estor)-log(lowbdor))/(sqrt(1/(sn011+0.5)+1/(sn101+0.5)+
                                          1/(sn012+0.5)+1/(sn102+0.5))/2);
    if ztest ge 1.96 then sumrej+1;
    end;  ** end of do nsimul;
   spower=sumrej/nsimul;
   aprxpower=spower;
put or 1-10 2 p011 11-20 2 p101 21-30 2 aprxn 31-35 aprxerr1 40-45 3 
    aprxpower 50-55 3;
end;  ** end of do p102 p012;
end;  ** end of do p101;
end;  ** end of do p011;
end;  ** end of do or;
end;  ** end of do lowbdor;