data step1;
** You may change the following number of simulations from 1000000 to 10000;
** to save the cpu time Based on our simulations, we find that the difference in;
** estimate of the final minimum number of patients under the same configuration;
** between these is within one or two cases at most;
nsimul=1000000;
** You may change the following alpha level and za to
the value as you wish;
alpha=0.025; za=1.96; zb=probit(0.8);
** When you change the following range for pxyz;
** you will need to assure that all the cell probabilities;
** fall 0 and 1 Otherwise, there will be 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 between two treatments;
do p011=0.05,0.10,0.20; ** the pi01 in group 1;
do p101=0.15,0.25; ** the po10 in group 1;
pc1=1-p011-p101; ** probability of concordance in group 1;
do p012=0.10; ** the pi01 in group 2;
p102=or**2*p101*p012/p011; ** the pi10 in group 2;
pc2=1-p012-p102; ** probability of concordance in group 2;
vlor=(1/p011+1/p101+1/p012+1/p102)/4; ** variance of log or estimate;
aprxn=int((za+zb)**2*vlor/(log(lowbdor)-log(or))**2)+1; ** approximated
sample size;
n=aprxn;
sumrej=0;
sumrej1=0;
do i=1 to nsimul;
sn011=0; sn101=0; sn012=0; sn102=0;
do nx=1 to n;
p102=lowbdor**2*p101*p012/p011; ** the pi10 in group 2;
pc2=1-p012-p102; ** probability of concordance in group 2;
x=rantbl(1432151,p011,p101,pc1); ** simulate the trinomial distribution
for group 1;
y=rantbl(1431171,p012,p102,pc2); *8 simulate the trinomial distribution
for group 2;
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;
sxplus=sn011+sn012; ** xplus;
sndis1=sn011+sn101; ** number of patients with discordance for group 1;
sndis2=sn012+sn102; ** number of patients with discordance for group 2;
if sxplus ge 1 and sndis1 ge 1 and sndis2 ge 1 then do;
a=max(0,sxplus-sndis2);
b=min(sndis1, sxplus);
numpowerx=0;
denpowerx=0;
nflag1=0;
nflag2=0;
do x1=a to b;
f12=lgamma(sndis1+1)-lgamma(x1+1)-lgamma(sndis1-x1+1)+
lgamma(sndis2+1)-lgamma(sxplus-x1+1)-lgamma(sndis2-sxplus+x1+1)+
x1*log(lowbdor**2);
nflag2=int(f12/50);
if nflag2 le nflag1 then do;
f12=f12-nflag1*50;
end;
if nflag2 gt nflag1 then do;
f12=f12-nflag2*50;
numpowerx=numpowerx*exp((nflag2-nflag1)*(-50));
denpowerx=denpowerx*exp((nflag2-nflag1)*(-50));
nflag1=nflag2;
end;
if x1 ge sn011 then do;
numpowerx+exp(f12); end;
denpowerx+exp(f12);
end; ** end of do x1=a to b;
powerx=numpowerx/denpowerx;
if powerx le alpha then sumrej1+1;
end; ** if sxplus ge 1 and sndis1 ge 1 and;
end; ** end of do i=1 to nsimul;
type1=sumrej1/nsimul;
aprxerr1=type1;
p102=or**2*p101*p012/p011; ** the pi10 in group 2;
pc2=1-p012-p102; ** probability of concordance in group 2;
do i=1 to nsimul;
sn011=0; sn101=0; sn012=0; sn102=0;
do nx=1 to n;
x=rantbl(1432151,p011,p101,pc1); ** simulate the trinomial distribution
for group 1;
y=rantbl(1431171,p012,p102,pc2); *8 simulate the trinomial distribution
for group 2;
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;
sxplus=sn011+sn012; ** xplus;
sndis1=sn011+sn101; ** number of patients with discordance for group 1;
sndis2=sn012+sn102; ** number of patients with discordance for group 2;
if sxplus ge 1 and sndis1 ge 1 and sndis2 ge 1 then do;
a=max(0,sxplus-sndis2);
b=min(sndis1, sxplus);
numpowerx=0;
denpowerx=0;
nflag1=0;
nflag2=0;
do x1=a to b;
f12=lgamma(sndis1+1)-lgamma(x1+1)-lgamma(sndis1-x1+1)+
lgamma(sndis2+1)-lgamma(sxplus-x1+1)-lgamma(sndis2-sxplus+x1+1)+
x1*log(lowbdor**2);
nflag2=int(f12/50);
if nflag2 le nflag1 then do;
f12=f12-nflag1*50;
end;
if nflag2 gt nflag1 then do;
f12=f12-nflag2*50;
numpowerx=numpowerx*exp((nflag2-nflag1)*(-50));
denpowerx=denpowerx*exp((nflag2-nflag1)*(-50));
nflag1=nflag2;
end;
if x1 ge sn011 then do;
numpowerx+exp(f12); end;
denpowerx+exp(f12);
end; ** end of do x1=a to b;
powerx=numpowerx/denpowerx;
if powerx le alpha then sumrej+1;
end; ** if sxplus ge 1 and sndis1 ge 1 and;
end; ** end of do i=1 to nsimul;
spower=sumrej/nsimul;
aprxpower=spower;
do while (spower lt 0.80);
n+5;
sumrej=0;
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;
sxplus=sn011+sn012;
sndis1=sn011+sn101;
sndis2=sn012+sn102;
if sxplus ge 1 and sndis1 ge 1 and sndis2 ge 1 then do;
a=max(0,sxplus-sndis2);
b=min(sndis1, sxplus);
numpowerx=0;
denpowerx=0;
nflag1=0;
nflag2=0;
do x1=a to b;
f12=lgamma(sndis1+1)-lgamma(x1+1)-lgamma(sndis1-x1+1)+
lgamma(sndis2+1)-lgamma(sxplus-x1+1)-lgamma(sndis2-sxplus+x1+1)+
x1*log(lowbdor**2);
nflag2=int(f12/50);
if nflag2 le nflag1 then do;
f12=f12-nflag1*50;
end;
if nflag2 gt nflag1 then do;
f12=f12-nflag2*50;
numpowerx=numpowerx*exp((nflag2-nflag1)*(-50));
denpowerx=denpowerx*exp((nflag2-nflag1)*(-50));
nflag1=nflag2;
end;
if x1 ge sn011 then do;
numpowerx+exp(f12); end;
denpowerx+exp(f12);
end; ** end of do x1=a to b;
powerx=numpowerx/denpowerx;
if powerx le alpha then sumrej+1;
end; ** if sxplus ge 1 and sndis1 ge 1 and;
end; ** end of do i=1 to nsimul;
spower=sumrej/nsimul;
end; ** end of do while (0.80 lt 0.80);
fpower=spower;
do while (spower gt 0.80);
n+(-1);
sumrej=0;
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;
sxplus=sn011+sn012;
sndis1=sn011+sn101;
sndis2=sn012+sn102;
if sxplus ge 1 and sndis1 ge 1 and sndis2 ge 1 then do;
a=max(0,sxplus-sndis2);
b=min(sndis1, sxplus);
numpowerx=0;
denpowerx=0;
nflag1=0;
nflag2=0;
do x1=a to b;
f12=lgamma(sndis1+1)-lgamma(x1+1)-lgamma(sndis1-x1+1)+
lgamma(sndis2+1)-lgamma(sxplus-x1+1)-lgamma(sndis2-sxplus+x1+1)+
x1*log(lowbdor**2);
nflag2=int(f12/50);
if nflag2 le nflag1 then do;
f12=f12-nflag1*50;
end;
if nflag2 gt nflag1 then do;
f12=f12-nflag2*50;
numpowerx=numpowerx*exp((nflag2-nflag1)*(-50));
denpowerx=denpowerx*exp((nflag2-nflag1)*(-50));
nflag1=nflag2;
end;
if x1 ge sn011 then do;
numpowerx+exp(f12); end;
denpowerx+exp(f12);
end; ** end of do x1=a to b;
powerx=numpowerx/denpowerx;
if powerx le alpha then sumrej+1;
end; ** if sxplus ge 1 and sndis1 ge 1 and;
end; ** end of do i=1 to nsimul;
spower=sumrej/nsimul;
if spower ge 0.80 then fpower=spower;
end; **do while power gt 0.80;
fin=n+1;
n=fin;
sumrej1=0;
do i=1 to nsimul;
sn011=0; sn101=0; sn012=0; sn102=0;
do nx=1 to n;
p102=lowbdor**2*p101*p012/p011; ** the pi10 in group 2;
pc2=1-p012-p102; ** probability of concordance in group 2;
x=rantbl(1432151,p011,p101,pc1); ** simulate the trinomial distribution
for group 1;
y=rantbl(1431171,p012,p102,pc2); *8 simulate the trinomial distribution
for group 2;
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;
sxplus=sn011+sn012; ** xplus;
sndis1=sn011+sn101; ** number of patients with discordance for group 1;
sndis2=sn012+sn102; ** number of patients with discordance for group 2;
if sxplus ge 1 and sndis1 ge 1 and sndis2 ge 1 then do;
a=max(0,sxplus-sndis2);
b=min(sndis1, sxplus);
numpowerx=0;
denpowerx=0;
nflag1=0;
nflag2=0;
do x1=a to b;
f12=lgamma(sndis1+1)-lgamma(x1+1)-lgamma(sndis1-x1+1)+
lgamma(sndis2+1)-lgamma(sxplus-x1+1)-lgamma(sndis2-sxplus+x1+1)+
x1*log(lowbdor**2);
nflag2=int(f12/50);
if nflag2 le nflag1 then do;
f12=f12-nflag1*50;
end;
if nflag2 gt nflag1 then do;
f12=f12-nflag2*50;
numpowerx=numpowerx*exp((nflag2-nflag1)*(-50));
denpowerx=denpowerx*exp((nflag2-nflag1)*(-50));
nflag1=nflag2;
end;
if x1 ge sn011 then do;
numpowerx+exp(f12); end;
denpowerx+exp(f12);
end; ** end of do x1=a to b;
powerx=numpowerx/denpowerx;
if powerx le alpha then sumrej1+1;
end; ** if sxplus ge 1 and sndis1 ge 1 and;
end; ** end of do i=1 to nsimul;
type1=sumrej1/nsimul;
ferr1=type1;
put or 1-10 2 p011 11-17 2 p101 18-24 2 aprxn 25-30 @33 "("
aprxerr1 34-38 3 @39 ")" @42 "[" aprxpower 43-47 3 @48 "]" fin 52-57 @60 "("
ferr1 61-65 3 @66 ")" @69 "[" fpower 70-74 3 @75 "]";
end; ** end of do p012;
end; ** end of do p101;
end; ** end of do p011;
end; ** end of do or;
end; ** end of do lowbdor;