Statistical power with SAS

/***************** powtest.sas *********************
   Non-centrality parameter must be given directly
   In this example, say mu1=12.5, mu2=13, mu3=18, mu4=21 
   and sigma = 2.5
****************************************************/
options linesize = 79 pagesize = 200;
data power;         /* Replace alpha, s, & p below        */
     alpha = 0.05;  /* Significance level                 */
     s = 3;         /* Numerator df = # IVs being tested  */
     p = 4;         /* There are p beta parameters        */
     n = 10  ;      /* Sample size                        */
     ncp = 18.16;     df2 = n-p;    oneminus = 1 - alpha;
     power = 1-probf(finv(oneminus,s,df2),s,df2,ncp);

     put '   ';
     put ' ***************************************************';
     put '   ';
     put '   For a multiple regression model with ' p 'betas, ';
     put '   testing ' s 'independent variables using ';
     put '   alpha = ' alpha ' and a sample size of ' n ',';
     put '   power = ' power ' for a non-centrality parameter';
     put '   value of phi = ' ncp;
     put '   ';
     put ' ****************************************************';
     put '   ';

Get output on the log file:

 ***************************************************

   For a multiple regression model with 4 betas,
   testing 3 independent variables using
   alpha = 0.05  and a sample size of 10 ,
   power = 0.7310757243  for a non-centrality parameter
   value of phi = 18.16

 ****************************************************
 

Here's a program for finding n to achieve desired power

 
 
/*****************  power.sas *********************/
options linesize = 79 pagesize = 200;
data fpower;        /* Replace alpha, s, p, and wantpow below    */
     alpha = 0.05;  /* Significance level                        */
     s = 6;         /* Numerator df = # IVs being tested         */
     p = 26;        /* There are p beta parameters               */
     ncp=18 ;       /* Non-centrality parameter phi              */
     wantpow = .90; /* Find n to yield this power.               */
     power = 0; n = p+1; oneminus = 1-alpha; /* Initializing ... */
     do until (power >= wantpow);
        df2 = n-p;
        power = 1-probf(finv(oneminus,s,df2),s,df2,ncp);
        n = n+1 ;
     end;
     n = n-1;
     put ' *********************************************************';
     put '   ';
     put '   For a multiple regression model with ' p 'betas, ';
     put '   testing ' s 'independent variables using alpha = ' alpha ',';
     put '   a sample size of ' n 'is needed';
     put '   in order to have probability ' wantpow 'of rejecting H0';
     put '   for a non-centrality parameter of phi = ' ncp ;
     put '   ';
     put ' *********************************************************';

Get output on the log file:

 *********************************************************

   For a multiple regression model with 26 betas,
   testing 6 independent variables using alpha = 0.05 ,
   a sample size of 218 is needed
   in order to have probability 0.9 of rejecting H0
   for a non-centrality parameter of phi = 18

 *********************************************************