data step1;
* This program is for calculation of the minimum required sample size;
* for assessing equivalence based on proportion ratio under a RCT;
* with noncompliance and missing outcomes;
* This program also calculates the simulated power for the corresponding;
* sample size so that we can evaluate the accuracy of the sample size estimate;
array spex1(j) sxpe1-sxpe2;
array spex2(j) sxpe3-sxpe4;
array spex3(j) sxpe5-sxpe6;
array spex(i) spex1-spex3;
array spey1(j) sype1-sype2;
array spey2(j) sype3-sype4;
array spey3(j) sype5-sype6;
array spey(i) spey1-spey3;
array sam(m) samx1-samx5;
array pow(m) epow1-epow5;
array fail(m) efail1-efail5;
za=1.645; ** the upper 5 percentile of the standard normal distribution;
nsimul=10000;
do rconc=0.30,0.70; ** conditional prob of response among compliers;
rcona=0.35; ** condtional prob of response for always-takers is;
** rcona*r considered here;
rconn=0.25; **conditional prob of response for never-takers;
put rconc 1-10 3;
do pattern= 1, 2;
put pattern 10-30;
do marprob=1 to 3;
if marprob eq 1 then do;
pc=0.95; pa=0.03; pn=0.02; ** proportion of subpopulations in different;
end; ** compliance categories;
if marprob eq 2 then do;
pc=0.80; pa=0.15; pn=0.05;
end;
if marprob eq 3 then do;
pc=0.60; pa=0.25; pn=0.15;
end;
delta=0.20; ** maximum acceptance level for equivalence;
desired=0.80; ** desired power;
do rr=0.90,1,1.1; *true ratio of prop of response between two treatments;
do numk=0.5,1,2,4,8; *ratio of sample allocation between two treatments;
if numk=0.5 then m=1;
if numk=1 then m=2;
if numk=2 then m=3;
if numk=4 then m=4;
if numk=8 then m=5;
if pattern eq 1 then do; ** prob of missing outcomes;
nomissc1=0.95; **for complier assigned to expermental treatment;
nomissc0=0.85; **for complier assigned to standard treatment;;
nomissa=0.90; **for always takers;
nomissn=0.80; **for never takers;
end;
if pattern eq 2 then do;
nomissc1=0.85;
nomissc0=0.75;
nomissa=0.80;
nomissn=0.70;
end;
xp1=pc*rconc*rr*nomissc1;
xp2=pa*rcona*rr*nomissa;
xp3=pn*rconn*nomissn;
xp4=pc*(1-rconc*rr)*nomissc1;
xp5=pa*(1-rcona*rr)*nomissa;
xp6=pn*(1-rconn)*nomissn;
xp7=pc*(1-nomissc1);
xp8=pa*(1-nomissa);
xp9=pn*(1-nomissn);
xpe1=xp1+xp2;
xpe2=xp3;
xpe3=xp4+xp5;
xpe4=xp6;
xpe5=xp7+xp8;
xpe6=xp9;
yp1=pc*rconc*nomissc0;
yp2=pa*rcona*rr*nomissa;
yp3=pn*rconn*nomissn;
yp4=pc*(1-rconc)*nomissc0;
yp5=pa*(1-rcona*rr)*nomissa;
yp6=pn*(1-rconn)*nomissn;
yp7=pc*(1-nomissc0);
yp8=pa*(1-nomissa);
yp9=pn*(1-nomissn);
ype1=yp1+yp3;
ype2=yp2;
ype3=yp4+yp6;
ype4=yp5;
ype5=yp7+yp9;
ype6=yp8;
n=5;
n1=n;
n0=int(n*numk+0.5);
v1=((xpe2+xpe4)*(1-(xpe2+xpe4))/n1+
(ype1+ype3)*(1-(ype1+ype3))/n0)/((ype1+ype3)-(xpe2+xpe4))**2;
v2=(xpe1*(1-xpe1)/n1+
ype2*(1-ype2)/n0)/(xpe1-ype2)**2;
v3=((xpe1+xpe3)*(1-(xpe1+xpe3))/n1+
(ype2+ype4)*(1-(ype2+ype4))/n0)/((ype2+ype4)-(xpe1+xpe3))**2;
v4=(xpe2*(1-xpe2)/n1+
ype1*(1-ype1)/n0)/(xpe2-ype1)**2;
c12=(ype2*(ype1+ype3)/n0+xpe1*(xpe2+xpe4)/n1)/
((ype1+ype3-(xpe2+xpe4))*(xpe1-ype2));
c13=((ype2+ype4)*(ype1+ype3)/n0+(xpe1+xpe3)*(xpe2+xpe4)/n1)/
((ype1+ype3-(xpe2+xpe4))*(xpe1+xpe3-(ype2+ype4)));
c14=((ype1-ype1*(ype1+ype3))/n0+(xpe2-xpe2*(xpe2+xpe4))/n1)/
((ype1+ype3-(xpe2+xpe4))*(ype1-xpe2));
c23=((ype2-ype2*(ype2+ype4))/n0+(xpe1-xpe1*(xpe1+xpe3))/n1)/
((xpe1-ype2)*(xpe1+xpe3-(ype2+ype4)));
c24=(ype1*ype2/n0+xpe1*xpe2/n1)/
((xpe1-ype2)*(ype1-xpe2));
c34=((ype2+ype4)*ype1/n0+(xpe1+xpe3)*xpe2/n1)/
((xpe1+xpe3-(ype2+ype4))*(ype1-xpe2));
var=v1+v2+v3+v4+2*c12-2*c13-2*c14-2*c23-2*c24+2*c34;
low=(log(1-delta)-log(rr))/sqrt(var)+za;
up=(log(1+delta)-log(rr))/sqrt(var)-za;
power=probnorm(up)-probnorm(low);
do while (power lt desired);
n+10;
n1=n;
n0=int(n*numk+0.5);
v1=((xpe2+xpe4)*(1-(xpe2+xpe4))/n1+
(ype1+ype3)*(1-(ype1+ype3))/n0)/((ype1+ype3)-(xpe2+xpe4))**2;
v2=(xpe1*(1-xpe1)/n1+
ype2*(1-ype2)/n0)/(xpe1-ype2)**2;
v3=((xpe1+xpe3)*(1-(xpe1+xpe3))/n1+
(ype2+ype4)*(1-(ype2+ype4))/n0)/((ype2+ype4)-(xpe1+xpe3))**2;
v4=(xpe2*(1-xpe2)/n1+
ype1*(1-ype1)/n0)/(xpe2-ype1)**2;
c12=(ype2*(ype1+ype3)/n0+xpe1*(xpe2+xpe4)/n1)/
((ype1+ype3-(xpe2+xpe4))*(xpe1-ype2));
c13=((ype2+ype4)*(ype1+ype3)/n0+(xpe1+xpe3)*(xpe2+xpe4)/n1)/
((ype1+ype3-(xpe2+xpe4))*(xpe1+xpe3-(ype2+ype4)));
c14=((ype1-ype1*(ype1+ype3))/n0+(xpe2-xpe2*(xpe2+xpe4))/n1)/
((ype1+ype3-(xpe2+xpe4))*(ype1-xpe2));
c23=((ype2-ype2*(ype2+ype4))/n0+(xpe1-xpe1*(xpe1+xpe3))/n1)/
((xpe1-ype2)*(xpe1+xpe3-(ype2+ype4)));
c24=(ype1*ype2/n0+xpe1*xpe2/n1)/
((xpe1-ype2)*(ype1-xpe2));
c34=((ype2+ype4)*ype1/n0+(xpe1+xpe3)*xpe2/n1)/
((xpe1+xpe3-(ype2+ype4))*(ype1-xpe2));
var=v1+v2+v3+v4+2*c12-2*c13-2*c14-2*c23-2*c24+2*c34;
low=(log(1-delta)-log(rr))/sqrt(var)+za;
up=(log(1+delta)-log(rr))/sqrt(var)-za;
power=probnorm(up)-probnorm(low);
end; ** end of do while;
btn=n-9;
ftn=n;
do n=btn to ftn by 1;
n1=n;
n0=int(n*numk+0.5);
v1=((xpe2+xpe4)*(1-(xpe2+xpe4))/n1+
(ype1+ype3)*(1-(ype1+ype3))/n0)/((ype1+ype3)-(xpe2+xpe4))**2;
v2=(xpe1*(1-xpe1)/n1+
ype2*(1-ype2)/n0)/(xpe1-ype2)**2;
v3=((xpe1+xpe3)*(1-(xpe1+xpe3))/n1+
(ype2+ype4)*(1-(ype2+ype4))/n0)/((ype2+ype4)-(xpe1+xpe3))**2;
v4=(xpe2*(1-xpe2)/n1+
ype1*(1-ype1)/n0)/(xpe2-ype1)**2;
c12=(ype2*(ype1+ype3)/n0+xpe1*(xpe2+xpe4)/n1)/
((ype1+ype3-(xpe2+xpe4))*(xpe1-ype2));
c13=((ype2+ype4)*(ype1+ype3)/n0+(xpe1+xpe3)*(xpe2+xpe4)/n1)/
((ype1+ype3-(xpe2+xpe4))*(xpe1+xpe3-(ype2+ype4)));
c14=((ype1-ype1*(ype1+ype3))/n0+(xpe2-xpe2*(xpe2+xpe4))/n1)/
((ype1+ype3-(xpe2+xpe4))*(ype1-xpe2));
c23=((ype2-ype2*(ype2+ype4))/n0+(xpe1-xpe1*(xpe1+xpe3))/n1)/
((xpe1-ype2)*(xpe1+xpe3-(ype2+ype4)));
c24=(ype1*ype2/n0+xpe1*xpe2/n1)/
((xpe1-ype2)*(ype1-xpe2));
c34=((ype2+ype4)*ype1/n0+(xpe1+xpe3)*xpe2/n1)/
((xpe1+xpe3-(ype2+ype4))*(ype1-xpe2));
var=v1+v2+v3+v4+2*c12-2*c13-2*c14-2*c23-2*c24+2*c34;
low=(log(1-delta)-log(rr))/sqrt(var)+za;
up=(log(1+delta)-log(rr))/sqrt(var)-za;
power=probnorm(up)-probnorm(low);
if round(power, 0.0000000001) ge desired then do;
fn=n;
n=ftn+100; end;
end; ** end of do n = btn to ftn;
n1=fn;
n0=int(n1*numk+0.5);
sumrej=0;
sumfail=0;
do numsim=1 to nsimul;
flag=0;
do i =1 to 3;
do j = 1 to 2;
spex=0; spey=0;
end;
end;
do numexp=1 to n1;
pcat=rantbl(1342151,xpe1,xpe2,xpe3,xpe4,xpe5,xpe6);
if pcat eq 1 then sxpe1+1/n1;
if pcat eq 2 then sxpe2+1/n1;
if pcat eq 3 then sxpe3+1/n1;
if pcat eq 4 then sxpe4+1/n1;
if pcat eq 5 then sxpe5+1/n1;
if pcat eq 6 then sxpe6+1/n1;
end;
do numexp=1 to n0;
pcat=rantbl(1342151,ype1,ype2,ype3,ype4,ype5,ype6);
if pcat eq 1 then sype1+1/n0;
if pcat eq 2 then sype2+1/n0;
if pcat eq 3 then sype3+1/n0;
if pcat eq 4 then sype4+1/n0;
if pcat eq 5 then sype5+1/n0;
if pcat eq 6 then sype6+1/n0;
end;
if round(sxpe1-sype2,0.00000000001) le 0 or
round((sype1+sype3)-(sxpe2+sxpe4),0.0000000001) le 0 or
round((sype1-sxpe2),0.000000001) le 0 or
round((sxpe1+sxpe3)-(sype2+sype4),0.0000000001) le 0 then
flag=1;
if round(sxpe1-sype2,0.00000000001) gt 0 and
round((sype1+sype3)-(sxpe2+sxpe4),0.0000000001) gt 0 and
round((sype1-sxpe2),0.000000001) gt 0 and
round((sxpe1+sxpe3)-(sype2+sype4),0.0000000001) gt 0 then do;
num=(sxpe1-sype2)*((sype1+sype3)-(sxpe2+sxpe4));
den=(sype1-sxpe2)*((sxpe1+sxpe3)-(sype2+sype4));
elrr=log(num/den);
v1=((sxpe2+sxpe4)*(1-(sxpe2+sxpe4))/n1+
(sype1+sype3)*(1-(sype1+sype3))/n0)/((sype1+sype3)-(sxpe2+sxpe4))**2;
v2=(sxpe1*(1-sxpe1)/n1+
sype2*(1-sype2)/n0)/(sxpe1-sype2)**2;
v3=((sxpe1+sxpe3)*(1-(sxpe1+sxpe3))/n1+
(sype2+sype4)*(1-(sype2+sype4))/n0)/((sype2+sype4)-(sxpe1+sxpe3))**2;
v4=(sxpe2*(1-sxpe2)/n1+
sype1*(1-sype1)/n0)/(sxpe2-sype1)**2;
c12=(sype2*(sype1+sype3)/n0+sxpe1*(sxpe2+sxpe4)/n1)/
((sype1+sype3-(sxpe2+sxpe4))*(sxpe1-sype2));
c13=((sype2+sype4)*(sype1+sype3)/n0+(sxpe1+sxpe3)*(sxpe2+sxpe4)/n1)/
((sype1+sype3-(sxpe2+sxpe4))*(sxpe1+sxpe3-(sype2+sype4)));
c14=((sype1-sype1*(sype1+sype3))/n0+(sxpe2-sxpe2*(sxpe2+sxpe4))/n1)/
((sype1+sype3-(sxpe2+sxpe4))*(sype1-sxpe2));
c23=((sype2-sype2*(sype2+sype4))/n0+(sxpe1-sxpe1*(sxpe1+sxpe3))/n1)/
((sxpe1-sype2)*(sxpe1+sxpe3-(sype2+sype4)));
c24=(sype1*sype2/n0+sxpe1*sxpe2/n1)/
((sxpe1-sype2)*(sype1-sxpe2));
c34=((sype2+sype4)*sype1/n0+(sxpe1+sxpe3)*sxpe2/n1)/
((sxpe1+sxpe3-(sype2+sype4))*(sype1-sxpe2));
var=v1+v2+v3+v4+2*c12-2*c13-2*c14-2*c23-2*c24+2*c34;
if (elrr-log(1.2))/sqrt(var) le -1*za and
(elrr-log(0.8))/sqrt(var) ge za then sumrej+1;
end; ** end of if gt 0 for all;
sumfail+flag;
end; ** end of do simul;
epower=sumrej/(nsimul-sumfail);
sam=n1;
pow=epower;
fail=sumfail/nsimul;
if m eq 5 then do;
put pc 1-5 2 pa 6-10 2 pn 11-15 2 rr 16-20 2
samx1 26-32 samx2 36-42 samx3 46-52 samx4 56-62 samx5 66-72/
@27 "(" epow1 28-32 3 @33 ")" @37 "(" epow2 38-42 3 @43 ")"
@47 "(" epow3 48-52 3 @53 ")" @57 "(" epow4 58-62 3 @63 ")"
@67 "(" epow5 68-72 3 @73 ")"/
@27 "[" efail1 28-32 3 @33 "]" @37 "[" efail2 38-42 3 @43 "]"
@47 "[" efail3 48-52 3 @53 "]" @57 "[" efail4 58-62 3 @63 "]"
@67 "[" efail5 68-72 3 @73 "]";
end; ** end of do if m eq 5;
end; ** end of do numk;
end; ** end of do rr;
end; ** end of do marprob;
end; ** do pattern;
end; ** end of do rconc;