data step1;
** This program calculates the approximate sample size and;
** its corresponding probability of correct selection based on;
** the product of the exact beta-binomial distribution of beta;
** You may change the parameter values to fit the situations you consider;

dd=0.05; ** [-dd, dd] is the ambiguous interval;
** You may replace 0.05 by any appropriate value as you wish;

za=1.282; ** This corresponds to 90% of the probability of correct selection;
** You may set za=1.645 for 95% of this probability;

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 'Approx' @10 'Cluster' @20 'Intra'
@30 'Response' @40 'Difference' @55 'Exact'/
@1 'sample' @10 'size' @20'corr' @30 'prob of' @40 'in prob' @55 'prob'/
@1 'size' @30 'treatment' @40 'between' @55 'of'/ @35 'I'
@40 'treatments' @55 'correct'/@55 'selection';

do pi1=0.20,0.30,0.40,0.50; ** This corresponds to the probability;
** of response in treatment 1;

delta=0.15; ** This denotes the true difference in the;
** probabilites of response between treatments;
pi2=pi1-delta;

do rho=0.10,0.20,0.50,0.80; ** This denotes the intraclass correlation;
t=(1-rho)/rho;
a1=t*pi1; b1=t*(1-pi1);
a2=t*pi2; b2=t*(1-pi2);

do r=1,2,4;
n=int(za**2*(pi1*(1-pi1)+pi2*(1-pi2))*(1+(r-1)*rho)/(r*(dd-delta)**2)+0.5);
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;
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;
end; ** end of do l=1;
end; ** end of do j=1;
put n 1-5 r 10-14 rho 19-23 2 pi1 30-36 2 delta 41-47 2 probdd 53-59 3;
end; ** end of do r=2;
end; ** end of do rho;
end; ** end of do pi1;