/******************************************************************************
 * Programme : Leistungsanalyse und Monte-Carlo-Simulation für Zwei-Mittelwert-Tests
 * Reference : LEISTUBACD
 * Source    : https://www.wearecas.eu/en/sampleCode/LEISTUBACD
 ******************************************************************************/

/* --- BLOC 1 --- */
proc power;
   twosamplemeans test=diff
      meandiff = 5 6
      stddev = 12 18
      alpha = 0.05 0.1
      ntotal = 100 200
      power = .;
run;

/* --- BLOC 2 --- */
ods graphics on;

proc power plotonly;
   twosamplemeans test=diff
      meandiff = 5 6
      stddev = 12 18
      alpha = 0.05 0.1
      ntotal = 100 200
      power = .;
   plot;
run;

ods graphics off;

   plot
      min=60
      yopts=(ref=0.9 crossref=yes)
      vary(color by stddev, linestyle by meandiff, symbol by alpha);

   ods output output=powdata;

/* --- BLOC 3 --- */
data tpow;
   do meandiff = 5, 6;
      do stddev = 12, 18;
         do alpha = 0.05, 0.1;
            do ntotal = 100, 200;
               ncp = ntotal * 0.5 * 0.5 * meandiff**2 / stddev**2;
               critval = finv(1-alpha, 1, ntotal-2, 0);
               power = sdf('f', critval, 1, ntotal-2, ncp);
               output;
            end;
         end;
      end;
   end;
run;

/* --- BLOC 4 --- */
proc print data=tpow;
run;

/* --- BLOC 5 --- */
%let meandiff =     5;
%let stddev   =    12;
%let alpha    =  0.05;
%let ntotal   =   100;
%let nsim     = 10000;

data simdata;
   call streaminit(123);
   do isim = 1 to ≁
      do i = 1 to floor(&ntotal/2);
         group = 1;
         y = rand('normal', 0        , &stddev);
         output;
         group = 2;
         y = rand('normal', &meandiff, &stddev);
         output;
      end;
   end;
run;

/* --- BLOC 6 --- */
ods exclude all;
proc ttest data=simdata;
   ods output ttests=tests;
   by isim;
   class group;
   var y;
run;
ods exclude none;

/* --- BLOC 7 --- */
data tests;
   set tests;
   where method="Pooled";
   issig = probt < α
run;

/* --- BLOC 8 --- */
proc freq data=tests;
   ods select binomial;
   tables issig / binomial(level='1');
run;