%macro power(version,
             data=,
              out=,
           effect=,
            calcs=,
               ss=ss3,
            alpha=0.05,
                n=,
            sigma=,
            delta=,
            debug=);

  %if &version ne %then %put POWER macro Version 1.2;

  %let opts = %sysfunc(getoption(notes)) validvarname=%sysfunc(getoption(validvarname));
  options nonotes 
         %if &sysver >= 7 %then %str(validvarname=V6);
  ;

  %let calcs=%upcase(&calcs);
  %let ss=%upcase(&ss);

* GET ERROR AND HYPOTHESIS DF AND
  SS FROM OUTSTAT DATA SET;
  data _null_;
    set &data(where=(_type_ in ("&ss",'ERROR')))
                end=lastobs;

    if _n_=1 then do;
      dfr=0;
      norig=0;
    end;

    norig+df;
    select(upcase(_source_));
      when ('ERROR') do;
        call symput('dfeorig',left(put(df,8.)));
        if ss gt 0 then
           call symput('sigorig',
              left(put(sqrt(ss/df),14.2)));
        else call symput('sigorig','');
      end;
      when (upcase("&effect")) do;
        dfr+df;
        call symput('dfh',left(put(df,8.)));
        call symput('ssh',left(put(ss,14.2)));
        call symput('fsamp',left(put(f,14.2)));
      end;
      otherwise do;
        dfr+df;
      end;
    end;

    if lastobs then do;
      call symput('dfr',left(put(dfr,8.)));
      call symput('norig',left(put((norig+1),8.)));
    end;
  run;

* PUT ORIGINAL DELTA INTO MACRO VARIABLE;
  data _null_;
    delta=sqrt(&ssh/&norig);
    call symput('delorig',left(put(delta,16.4)));
    stop;
  run;

* PREPEND VALUES FOR N SIGMA DELTA THAT
  OCCUR IN DATA TO USER-SPECIFIED VALUES;
  %let n=&norig &n;
  %let sigma=&sigorig σ
  %let delta=&delorig δ
  %if &debug ne %then %do;
   %put n=&n;
   %put sigma=σ
   %put delta=δ
   %put alpha=α
   %put dfeorig=&dfeorig;
   %put sigorig=&sigorig;
   %put dfr=𝔡
   %put dfh=&dfh;
   %put ssh=&ssh;
   %put fsamp=&fsamp;
  %end;


* MACRO TO ADD COMMAS AFTER VALUES
  SPECIFIED FOR ALPHA, N, SIGMA, DELTA;
  %macro comma(string,varname);
    %*global &varname;
    %local count word;
    %let count=1;
    %let word=%qscan(&string,&count,%str( ));
    %let &varname=&word;
    %do %while (&word ne);
       %let count=%eval(&count+1);
       %let word=%qscan(&string,&count,%str( ));
       %if (&word ne) %then
           %let &varname=&&&varname, &word;
    %end;
    %if &debug ne %then %put &varname = &&&varname;
  %mend comma;


* CREATE COMMA DELIMITTED LIST FOR USER-SPECIFIED
  VALUES OF ALPHA, N, SIGMA, DELTA;
  %let given=alpha n sigma delta;

  %do i=1 %to 4;
     %let current=%scan(&given,&i,%str( ));
     %comma(&&&current,&current);
  %end;

options notes;
* PERFORM CALCULATIONS;
  data &out;
    do alpha=α
      do number=&n;
         dfe=number-&dfr-1;
         do sigma=σ
            do delta=δ

          * CALCULATE LAMBDA AND POWER;
            lambda=(number*delta**2)/(sigma**2);
            astar=1-alpha;
            fcrit=finv(astar,&dfh,dfe);
            if lambda>135 then power=1.0;
            else power=
               1-probf(fcrit,&dfh,dfe,lambda);


          * CALCULATE ADJUSTED LAMBDA AND
            CI FOR ORIGINAL DELTA;
            if delta=&delorig and
               index("&calcs",'ADJPOW')>0 then do;
               lamadj=max(0,(lambda*(&dfeorig-2)/&dfeorig)-&dfh);
               if lamadj>135 then adjpow=1.0;
               else adjpow=1-probf(fcrit,&dfh,dfe,lamadj);
            end;
            else if delta=&delorig then adjpow=.U;
            else do;
               adjpow=.N;
            end;


          * GET CI ON LAMBDA AND POWER;
            if index("&calcs",'POWCI')>0 and
               adjpow not in (.N,.U) then do;
              %if &fsamp=I %then %do;
                powlow=.U; powup=.U;
              %end;
              %else %do;
               lamlow=&dfh*(max(0,(sqrt(&fsamp)-sqrt(fcrit))))**2;
               if lamlow>135 then powlow=1.0;
               else powlow=1-probf(fcrit,&dfh,dfe,lamlow);

               lamup=&dfh*(sqrt(&fsamp)+sqrt(fcrit))**2;
               if lamup>135 then powup=1.0;
               else powup=1-probf(fcrit,&dfh,dfe,lamup);
              %end;
            end;
            else do;
               powlow=.N;
               powup=.N;
            end;

          * FIND LEAST SIGNIFICANT N;
            if number=&norig and
               index("&calcs",'LSN')>0 then do;
               niter=&dfr+2;
               lstar=(niter*delta**2)/(&dfh*sigma**2);

               do until (diff<0.0000001);
                 niter=niter+1;
                 errn=niter-&dfr-1;
                 lstar=(niter*delta**2)/(&dfh*sigma**2);
                 diff=astar-probf(lstar,&dfh,errn);
               end;

               lsn=niter;
               lsndfe=lsn-&dfr-1;
               lamblsn=(lsn*delta**2)/(sigma**2);
               astar=1-alpha;
               fcrit=finv(astar,&dfh,lsndfe);
               if lamblsn>135 then powlsn=1.0;
               else powlsn=1-probf(fcrit,&dfh,lsndfe,lamblsn);
            end;
            else do;
               lsn=.N;
               powlsn=.N;
            end;

            output;

            end; *delta;
         end;    *sigma;
      end;       *number;
    end;         *alpha;

    label  alpha='Type I Error Rate'
          number='Sample Size'
           sigma='Root Mean Square Error'
           delta='Effect Size'
           power='Power of Test'
          adjpow='Adjusted Power'
          powlow='Confidence Interval: Lower Limit'
           powup='Confidence Interval: Upper Limit'
             lsn='Least Significant Number'
          powlsn='Power when N=LSN';

    keep alpha number sigma delta
       %if %index(&calcs,POWER)>0 %then %str(power);
       %if %index(&calcs,ADJPOW)>0 %then %str(adjpow);
       %if %index(&calcs,POWCI)>0 %then %str(powlow powup);
       %if %index(&calcs,LSN)>0 %then %str(lsn powlsn);;
  run;

options nonotes;
  %let ssnum=%substr(&ss,3,1);
* PRINT RESULTS;
  proc print data=&out noobs label;
    var alpha number sigma delta
        %if %index(&calcs,POWER)>0 %then %str(power);
        %if %index(&calcs,ADJPOW)>0 %then %str(adjpow);
        %if %index(&calcs,POWCI)>0 %then %str(powlow powup);
        %if %index(&calcs,LSN)>0 %then %str(lsn powlsn);;
    title1 "Power Calculation for effect %upcase(&effect)";
    title2 "Type &ssnum Sums of Squares";
  run;
  title;
options &opts;
%mend power;