/******************************  matpow2.sas ******************************/
title 'Power analysis for fixed-effects factorial ANOVA: Matrix approach';
title2 'Specify a beta vector instead of just an effect';
proc iml;
/********* Edit this input: Rows of matrices are separated by commas  ********/
     alpha = 0.05;
     wantpow = .80;
     f = {1,1,1,1,1,1};          /* Relative sample sizes       */
     C = { 1 -1 -1  1  0  0,     /* Contrast matrix             */
           0  0  1 -1 -1  1};
     beta = {0,.25,0,.25,0,-.25};

              /* In standard deviation units */
/*****************************************************************************/
     eff = C*beta;
     r = nrow(f) ; q = nrow(eff); f = f/sum(f);
/*** Echoing input ***/
     print " Significance level: " alpha;
     print " Desired power: " wantpow;
     print " Relative sample sizes:" f;
     print " Contrast matrix "; print C;
     print " Beta matrix (in SD units): " beta;
     print " Effect (h) matrix in SD units: " eff;
/*** Checking input ***/
     print " Checking input ...  ";
     inerr = 0; /* Input error = no. This may be changed below */
     if nrow(C) ^= q then do; 
        print "Error: Length of eff must equal number of rows in C";
        inerr=1;
     end;
     if ncol(C) ^= r then do;
        print "Error: Length of f must equal number of cols in C";
        inerr=1;
     end;
     if eff`*eff = 0 then do;
        print "Error: eff can't be zero; this would cause an infinite loop!";
        inerr=1;
     end;
     aa = min(f##2);
     if aa = 0 then do;
        print "Error: No sample size may be zero.";
        inerr=1;
     end;

/*** Proceed if no input errors ***/
     if inerr = 1 then abort;
        else do;
        print "Input appears to be okay.  ";
             core = inv(C*inv(diag(f))*C`);
             effsize = eff`*core*eff;  print " Effect size = " effsize;
             power = 0; n = r+1; oneminus = 1-alpha;  /* Initializing ...*/
             do until (power >= wantpow);
                n = n+1 ;
                ncp = n * effsize;
                df2 = n-r;
                power = 1-probf(finv(oneminus,q,df2),q,df2,ncp);
             end; /*  End Loop */
             print " Required sample size is " n;
        end; /* End computation (Conditional on no input errors) */
quit;