/* Program to use simulation to compute te power of the t test 
   for mean differences in (0 to 2 by 0.1) */
%macro ODSOff(); /* Call prior to BY-group processing */
   ods graphics off;
   ods exclude all;
   ods noresults;
%mend;
 
%macro ODSOn(); /* Call after BY-group processing */
   ods graphics on;
   ods exclude none;
   ods results;
%mend;
 
%let n1 = 10;           /* group sizes*/
%let n2 = 10;
%let NumSamples = 5000; /* number of simulated samples */  

data PowerSim(drop=i);
call streaminit(321);
do Delta = 0 to 2 by 0.1;
   do SampleID = 1 to &NumSamples;
      c = 1;                            /* Group 1 */
      do i = 1 to &n1;
         x = rand("Normal", 0, 1);      /* x ~ N(0,1) */
         output;
      end;
      c = 2;                            /* Group 2 */
      do i = 1 to &n2;
         x = rand("Normal", Delta, 1); /* x ~ N(Delta, 1) */
         output;
      end;
   end;
end;
run;

/* 2. Compute (pooled) t test for each sample */
%ODSOff 
proc ttest data=PowerSim; 
   by Delta SampleID; 
   class c; 
   var x; 
   ods output ttests=TTests(where=(method="Pooled")); 
run; 
%ODSOn
 
/* Construct indicator var for obs that reject H0 at 0.05 significance */ 
data Results; 
   set TTests; 
   RejectH0 = (Probt <= 0.05); 
run;

/* 3. Compute proportion: (# that reject H0)/NumSamples and CI */ 
%ODSoff
proc freq data=Results; 
   by Delta;
   tables RejectH0 / nocum binomial(level='1');
   output out=Est binomial;
run;
%ODSOn

/* for comparison, compute exact power */
proc power;
  twosamplemeans  power = .           /* missing ==> "compute this" */
    meandiff= 0 to 2 by 0.1           /* delta = 0, 0.1, ..., 2     */
    stddev=1                          /* N(delta, 1)                */
    ntotal=20;                        /* 20 obs in the two samples  */
  ods output Output=Power;            /* output results to data set */
run;

data Combine;
   set Est Power;
run;

/* overlay the simulated results with the cuve of exact power */
title "Power of the t Test";
title2 "Samples are N(0,1) and N(delta,1), n1=n2=10";
proc sgplot data=Combine noautolegend;
   series x=MeanDiff y=Power;
   scatter x=Delta y=_BIN_ / yerrorlower=L_Bin yerrorupper=U_Bin;
   inset ("Number of Samples"="&NumSamples") / border;
   yaxis min=0 max=1 label="Power (1 - P[Type II Error])" grid;
   xaxis label="Difference in Population Means" grid;
run;