Retrospective power analysis for anova

/*  retrospective power analysis for anova              */
/*  stata version 1.0 --6/23/09 -- pbe                  */
/*                                                      */
/* acknowledgement,  code adapted from                  */
/* rpower.sas by Michael Friendly  <friendly@yorku.ca>  */
/* version 1.1 november 2000                            */
/*                                                      */
/* instructions:                                        */
/* step 1 -- run anova model                            */
/* step 2 -- run test command on factor of interest     */
/* step 3 -- run code fragment (change alpha as needed) */
/*                                                      */
/* This do-file does not advocate retrospective power   */
/* analysis for anova but acknowledges the reality that */
/* retrospective power analysis is often demanded of    */
/* students and/or authors                              */

version 9.2
   local alpha = .05
   local sse=e(rss)
   local dfe=e(df_r)
   local mse=`sse'/`dfe'
   local ss=r(ss)
   local df =r(df)
/* calculate F non-centrality (Wright & O'Brien, 1988, p.1098) */
   local nc = `ss'/`mse'
/* adjusted F non-centrality, correcting for bias */
   local nca= max( 0, ( `nc'*(`dfe'-2)/`dfe')-`df' )
/* compute critical value of F */
   local fcrit = invFtail(`df', `dfe', `alpha')
/* find power from non-centrality */
   local power = nFtail(`df', `dfe', `nc', `fcrit')
   local powera= nFtail(`df', `dfe', `nca',`fcrit')
   display "  nominal power = " %5.4f `power'
   display " adjusted power = " %5.4f `powera'