data step1;

/* This program is for calculating Type I error of the proposed test procedures for testing

the homogeneity of risk difference under balanced design. In the program readers may modify

the assigned value in the “do” statement according to your interest. Note that since the probability

must fall between 0 and 1, the sum of assigned values for “con” and “dd” cannot be greater than 1. */

options ls=80;

array x(i) xx1-xx300;

array y(i) yy1-yy300;

array tx(j) txx1-txx300;

array ty(j) tyy1-tyy300;

array terrz(l) errz1-errz4;

array terrs(l) errs1-errs4;

perct=probnorm(1.645);

za=probit(0.95);

do con=0.20,0.50; ** the constant c defined for the range of p0k;

prob=0.50; ** the parameter pi for the beta distribution;

tab=2; ** the parameter T for the beta distribution;

** T=2 and pi=0.50 -- uniform distribution;

a=prob*tab; b=prob*tab;

do dd=0.0,0.10,0.20,0.40; ** the underlying common risk difference;

do n=30,50,100,200; ** the number of strata;

do m=2,3,5,10; ** the number of subjects per stratum in each treatment;

sumf1=0; sumf2=0;

fail1=0;

nsimul=1000; ** number of repeated samples for simulations;

bsimul=200; ** number of bootstrap samples;

sumest=0;

numer=0;

do nsi=1 to nsimul;

do j =1 to n;

tx=0; ty=0;

end;

sumden=0; sumnum=0;

sumx=0; sumy=0;

do j=1 to n;

aa=rangam(1342111,a); bb=rangam(1342111,b);

do while (round(aa,0.0000000000000001) le 0 or

round(bb,0.0000000000000001) le 0);

aa=rangam(1342111,a); bb=rangam(1342111,b);

end;

cc=aa/(aa+bb);

pi2=con*cc;

pi1=pi2+dd;

tx=ranbin(1342121,m,pi1);

ty=ranbin(1342121,m,pi2);

sumx+tx;

sumy+ty;

sumnum+(tx/m-ty/m);

sumden+1;

end;

ediff=sumnum/sumden;

num1=0;

do j = 1 to n;

num1+(tx/m-ty/m-ediff)**2/n;

end;

ex1=((sumx/(m*n))-((sumy/(n*m))**2/3+(sumx/(m*n))**2))/m;

ex2=((sumy/(n*m))-4*(sumy/(n*m))**2/3)/m;

num3=ex1+ex2;

numo=num1-num3;

nums=0;

numss=0;

sumgt=0;

do bi=1 to bsimul;

sumnum=0;

sumden=0;

sumx=0;

sumy=0;

do i=1 to n;

j=int(ranuni(1332111)*n)+1;

x=tx;

y=ty;

sumx+x;

sumy+y;

sumnum+(x/m-y/m);

sumden+1;

end; ** end of do i=1 to n;

ediff=sumnum/sumden;

num1=0;

do i=1 to n;

num1+(x/m-y/m-ediff)**2/n;

end;

ex1=((sumx/(m*n))-((sumx/(m*n))**2+(sumy/(n*m))**2/3))/m;

ex2=((sumy/(n*m))-4*(sumy/(n*m))**2/3)/m;

num3=ex1+ex2;

nums+(num1-num3);

if (num1-num3) ge 0 then sumgt+1;

numss+(num1-num3)**2;

end; ** end of do bi=1 to bsimul;

estvar=(numss-nums**2/bsimul)/(bsimul-1);

den=estvar;

if den gt 0 then do;

t=numo/sqrt(den);

if t gt za then sumf1+1;

end;

if den le 0 then fail1+1;

if sumgt ge (bsimul*perct) then sumf2+1;

end; ** end of do nsimul = 1 to 10000;

if m eq 2 then l=1;

if m eq 3 then l=2;

if m eq 5 then l=3;

if m eq 10 then l=4;

error1=sumf1/(nsimul-fail1);

error2=sumf2/nsimul;

terrz=error1;

terrs=error2;

profail=fail1/nsimul;

if l eq 4 then do;

put con 1-5 2 dd 6-10 2 n 11-15 errz1 16-21 3 errs1 22-27 3

errz2 28-33 3 errs2 34-39 3 errz3 40-45 3 errs3 46-51 3

errz4 52-57 3 errs4 58-63 3;

end; ** end of if l eq 4;

end; ** end of do m;

end; ** end of do n;

end; ** end of do dd;

end; ** end of do con;