/***********************************************************/
/* This macro calculates powers via Monte Carlo.  Arguments*/
/* are number of unexposed and exposed, and disease probab-*/
/* ilities for both groups.  Macro then makes a data set   */
/* with 40000 records.  Each group of 4 corresponds to a   */
/* sampled 2x2 table, with pairs of rows representing      */
/* independent binomial trials.                            */
/***********************************************************/
%macro dosamp(n1,n2,p1,p2);
   data rand; retain y seed; seed=0;
      do i=1 to 10000; 
      gp1=0; gp2=0; x1=ranbin(seed,&n1,&p1); y=x1; output;
      gp1=0; gp2=1; x1=&n1-y; output;
      gp1=1; gp2=0; x1=ranbin(seed,&n2,&p2); y=x1; output;
      gp1=1; gp2=1; x1=&n2-y; output;
      end; run;
   proc freq data=rand noprint; output out=outset chisq; 
      table gp1*gp2/chisq; weight x1; by i ; run;
   data outset; set outset; fisher2=0; fisher1=0; chisq=0; 
      if xpr_fish<.025 then fisher2=1; 
      if xp2_fish<.05 then fisher1=1; 
      if P_PCHI <.05 then chisq=1; run;
   proc means data=outset noprint; 
      output out=power mean=;
      var fisher2 fisher1 chisq; run;
   data power; set power; drop _freq_ _type_; 
      p0=(&p1+&p2)/2;
      se0=sqrt(p0*(1-p0)/&n1+p0*(1-p0)/&n2);
      seA=sqrt(&p1*(1-&p1)/&n1+&p2*(1-&p2)/&n2);
      mA=abs(&p2-&p1);
      aprx1=1-probnorm((1.96*se0-mA)/seA)
         /* +probnorm((-1.96*se0-mA)/seA) */;
      se0=sqrt(1/(4*&n1)+1/(4*&n2)); seA=se0;
      mA= abs(arsin(sqrt(&p2))-arsin(sqrt(&p1)));
      aprx2=1-probnorm((1.96*se0-mA)/seA)
         /* +probnorm((-1.96*se0-mA)/seA) */;
      n1=&n1; n2=&n2; p1=&p1; p2=&p2; drop mA p0 se0 seA;
      run;
   data powerout; set powerout power; run;
%mend dosamp;
data powerout; run;
%dosamp( 10, 10,.6,.4); %dosamp( 15,  5,.6,.4);
%dosamp( 10, 10,.8,.2); %dosamp( 20, 10,.6,.4);
%dosamp(200,100,.6,.4);
proc print noobs data=powerout; run;