data step1;
** This program calculates the minimum required number of clusters;
** for a given cluster size and a desired probability;
** of correct selection on the basis of the product of the exact
** beta-binomial distributions for the case of two arms;
dd=0.05; ** This denotes the ambiguous interval [-dd, dd];
lambda0=0.90; ** This denotes the desired probability of correction;
** selection;
array probt(i) prot1-prot6;
array probx(j) prox1-prox2000;
array proby(l) proy1-proy2000;
array sumx(k) sum1-sum2000;
array sumy(k) su1-su2000;
put @1 'prob of' @11 'prob diff' @21 'intra' @31 'cluster' @41 'required'
@51 'exact'/
@1 'response' @11 'between' @21 'correln' @31 'size' @41 'sample'
@51 'prob of'/ @1 'treat I' @11 'treatments' @41 'size' @51 'corr select';
do pi1=0.20 to 0.50 by 0.10; ** This represents the probability of response;
** for treatment 1;
** You may change this to fit the situation;
** which you wish;
delta=0.15; ** This denotes the difference between the;
** probability of response between treatment 1;
** and treatment 2;
pi2=pi1-delta;
do rho=0.10, 0.20,0.50,0.80; ** This is the intraclass correlation;
** You may change this to fit the situations;
** which you consider;
t=(1-rho)/rho;
a1=t*pi1; b1=t*(1-pi1);
a2=t*pi2; b2=t*(1-pi2);
do r=1,2,4; ** This is the cluster size;
** You may change this to fit the situation;
** which you consider;
n=10;
do k=1 to 2000; sumx=0; sumy=0; end;
do j=1 to 2000; probx=0; end;
do l=1 to 2000; proby=0; end;
do xi=0 to r;
i=xi+1;
probt=exp(lgamma(t)+lgamma(xi+a1)+lgamma(r+b1-xi)+lgamma(r+1)-
(lgamma(xi+1)+lgamma(r-xi+1)+lgamma(r+t)+lgamma(a1)+lgamma(b1)));
end;
do n1=1 to n;
if n1 = 1 then do;
do j=1 to r+1;
i=j;
probx=probt;
end;
end; ** end of if n=1;
if n1 ge 2 then do;
do i=1 to r+1;
do j=1 to (n1-1)*r+1;
k=(i-1)+j;
sumx+probx*probt;
end; ** end of do j=1 to (n1-1)*r+1;
end; ** end of do i=1 to r+1;
do j=1 to n1*r+1;
k=j;
probx=sumx;
sumx=0;
end;
end; ** end of if n1 ge 2;
end; ** end of do n1=1 to n;
do yi=0 to r;
i=yi+1;
probt=exp(lgamma(t)+lgamma(yi+a2)+lgamma(r+b2-yi)+lgamma(r+1)-
(lgamma(yi+1)+lgamma(r-yi+1)+lgamma(r+t)+lgamma(a2)+lgamma(b2)));
end;
do n2=1 to n;
if n2 = 1 then do;
do l=1 to r+1;
i=l;
proby=probt;
end;
end; ** end of if n2=1;
if n2 ge 2 then do;
do i=1 to r+1;
do l=1 to (n2-1)*r+1;
k=(i-1)+l;
sumy+proby*probt;
end; ** end of do l=1 to (n2-1)*r+1;
end; ** end of do i=1 to r+1;
do l=1 to n2*r+1;
k=l;
proby=sumy;
sumy=0;
end;
end; ** end of if n2 ge 2;
end; ** end of do n2=1 to n;
probdd=0;
probaa=0;
probda=0;
do j=1 to n*r+1;
do l=1 to n*r+1;
diff=(j-l)/(n*r);
if diff gt dd then probdd+probx*proby;
if diff gt dd then probda+probx*proby;
if -dd le diff le dd then probaa+probx*proby;
if -dd le diff le dd then probda+0.5*(probx*proby);
end; ** end of do l=1;
end; ** end of do j=1;
do while (probda lt lambda0); ** This incorporates the probability of;
** of correct selection due to the;
** ambiguous cases based on formula 8;
** do while (probdd lt lambda0); ** This does not include the ambiguous;
** cases when calculating the probability;
** of correct selection based on formula 5;
** If you want to apply this, you can remove;
** the two ** from this statement and;
** put one or two ** in the previous;
** do while statement;
** Also, don't forget changing the;
** output statement at the end of
** program as well;
n+1;
do k=1 to 2000; sumx=0; sumy=0; end;
do j=1 to 2000; probx=0; end;
do l=1 to 2000; proby=0; end;
do xi=0 to r;
i=xi+1;
probt=exp(lgamma(t)+lgamma(xi+a1)+lgamma(r+b1-xi)+lgamma(r+1)-
(lgamma(xi+1)+lgamma(r-xi+1)+lgamma(r+t)+lgamma(a1)+lgamma(b1)));
end;
do n1=1 to n;
if n1 = 1 then do;
do j=1 to r+1;
i=j;
probx=probt;
end;
end; ** end of if n=1;
if n1 ge 2 then do;
do i=1 to r+1;
do j=1 to (n1-1)*r+1;
k=(i-1)+j;
sumx+probx*probt;
end; ** end of do j=1 to (n1-1)*r+1;
end; ** end of do i=1 to r+1;
do j=1 to n1*r+1;
k=j;
probx=sumx;
sumx=0;
end;
end; ** end of if n1 ge 2;
end; ** end of do n1=1 to n;
do yi=0 to r;
i=yi+1;
probt=exp(lgamma(t)+lgamma(yi+a2)+lgamma(r+b2-yi)+lgamma(r+1)-
(lgamma(yi+1)+lgamma(r-yi+1)+lgamma(r+t)+lgamma(a2)+lgamma(b2)));
end;
do n2=1 to n;
if n2 = 1 then do;
do l=1 to r+1;
i=l;
proby=probt;
end;
end; ** end of if n2=1;
if n2 ge 2 then do;
do i=1 to r+1;
do l=1 to (n2-1)*r+1;
k=(i-1)+l;
sumy+proby*probt;
end; ** end of do l=1 to (n2-1)*r+1;
end; ** end of do i=1 to r+1;
do l=1 to n2*r+1;
k=l;
proby=sumy;
sumy=0;
end;
end; ** end of if n2 ge 2;
end; ** end of do n2=1 to n;
probdd=0;
probaa=0;
probda=0;
do j=1 to n*r+1;
do l=1 to n*r+1;
diff=(j-l)/(n*r);
if diff gt dd then probdd+probx*proby;
if diff gt dd then probda+probx*proby;
if -dd le diff le dd then probaa+probx*proby;
if -dd le diff le dd then probda+0.5*probx*proby;
end; ** end of do l=1;
end; ** end of do j=1;
end; ** end of do while;
put pi1 1-5 2 delta 11-15 2 rho 21-25 2 r 31-35 n 41-45 probda 51-56 3;
** This output statement is for the case in which we incorporate;
** the ambiguous cases into the probability of correction;
**put pi1 1-5 2 delta 11-15 2 rho 21-25 2 r 31-35 n 41-45 probdd 51-56 3;
** This output statement is for the case where we do not include;
** the ambiguous cases when calculate probability of correction;
** selection;
** When wishing this output statement, you can remove the two **;
** from this put statement and put one or two ** ahead of the previous;
** put statement to disable the previous one;
end; ** end of do r=2;
end; ** end of do rho;
end; ** end of do pi1;