/* sinfm5.sas                            */
 /*                                       */
 /* compute power of test to detect trend */
 /*  of infant mortality and income       */
 /*                                       */
  options ls=80 ;
 /* read county economic data             */
  data econ replicate ;
    length county $12 ;
    drop x1-x7 ;
    infile 'NCCounty03q4.txt' dlm='09'x firstobs=2 ;
    seed = 5151917 ;                       * set seed ;
 /*                              ***** skip some vars  */
    input fips county cpop cpopr x1-x7 hinc pcpi90 pcpi01 pov ;
 /* now just the bigger counties           */
 /*  but avoid the two biggest: Wake & Mecklenburg */
    if( cpopr < 3 ) then delete ;
    if( cpopr > 75 ) then delete ;
    output econ ;                          * old econ dataset ;
 /* generate some new y's                  */
    do slope = 0.000 to 0.0005 by 0.0001 ; * outer loop ;
    do rep = 1 to 1000 ;                    * 1000 reps ;
      totrate = 12.5 - slope*pcpi01 + 2*rannor(seed) ; * new y ;
      output replicate ;
    end ;
    end ;
 run ;
 /* should have 73*1000*6 obs              */
  proc sort data=replicate ;
    by slope rep ;                         * now two variables ;
 run ;
 /* take a look                            */
  proc print data=replicate (obs=5) ;
    title 'top of replicate dataset' ;
 run ;
  proc print data=replicate (firstobs=437995) ;
    title 'bottom of replicate dataset' ;
 run ;
 /* now regression analysis                */
 /*                                        */
 /*                                        */
 /*   ******* close   means turn off usual destination(listing) */
  ods listing close ;                      
  proc reg data=replicate ;                * no need for noprint ;
    model totrate = pcpi01 ;               * income ;
    ods output ParameterEstimates=pe ;     * put info into dataset pe ;
    title 'regression analysis of infant mortality on income' ;
    by slope rep ;                         * BY variable ;
 run ;
  ods listing ;
 /*   ******* go back to usual printing     */
   proc print data=pe (obs=10) ;
     title 'top of pe dataset created by ODS' ;
 run ;
   proc print data=pe (firstobs=1190) ;
     title 'bottom of pe dataset created by ODS' ;
 run ;
 /* analyze output dataset                  */
 /*                                         */
 /* simple way to count                     */
   data newpe ;
     set pe ;
     if( variable eq 'Intercept' ) then delete ;
     reject = (tvalue < -1.667 ) ;
 run ;
   proc means data=newpe mean nway noprint ;
      var reject ;
      class slope ;
      output out=b mean=power ;
      title 'fraction of replications rejecting' ;
 run ;
   proc print data=b ;
     title 'power dataset' ;
 run ;