/*********************************  matpow1.sas ******************************/
title 'Power analysis for fixed-effects factorial ANOVA: Matrix approach';
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};
     eff = {0,-.5};              /* In standard deviation units */
/*****************************************************************************/
     p = 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 " Effect 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) ^= p 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 = p; oneminus = 1-alpha;  /* Initializing ...*/
             do until (power >= wantpow);
                n = n+1 ;
                ncp = n * effsize;
                df2 = n-p;
                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;