Simultaneous Bonferroni Confidence Intervals

 /****************************************************************/
 /*       S A S   S A M P L E   L I B R A R Y                    */
 /*                                                              */
 /*    NAME: BONFER                                              */
 /*   TITLE: Simultaneous Bonferroni Confidence Intervals        */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Simultaneous Confidence Intervals,                  */
 /*   PROCS: MEANS PRINT SHEWHART                                */
 /*    DATA:                                                     */
 /*                                                              */
 /*     REF: R.G. Miller (1980), Simultaneous Statistical        */
 /*          Inference, New York: Springer-Verlag.               */
 /*                                                              */
 /*          Robert N. Rodriguez (1996), "Health Care            */
 /*          Applications of Statistical Process Control:        */
 /*          Examples Using the SAS System", Procedings of the   */
 /*          Twenty-First Annual SAS Users Group International   */
 /*          Conference, 1381-1396.                              */
 /*                                                              */
 /*   NOTES:                                                     */
 /*                                                              */
 /****************************************************************/

 options ls=64 nodate nonumber;
 goptions htext= 3.5 pct htitle=4.0 pct
          ftext='albany amt';

 options ls=64 nodate nonumber;

 * Simultaneous Bonferroni confidence intervals;

 data csection;
    length id $ 2 n94 n95 $4;
    input id csect94 vag94 csect95 vag95;
    total94 = csect94 + vag94;
    total95 = csect95 + vag95;
    n94     = left( put( total94, 4. ) );
    n95     = left( put( total95, 4. ) );
    label id  = 'Medical Group Identification Number'
          n94 = 'Number of Deliveries'
          n95 = 'Number of Deliveries';
    cards;
 1A      163     907     150     773
 1K       55     314      45     253
 1B       52     179      34     136
 1D       19     128      18     114
 3I       21      98      20      86
 3M       15      81      12      93
 1E       15      52      10      67
 1N        6      43      19      55
 1Q       12      67       7      62
 3H        7      65      11      54
 1R        4      43      11      38
 1H        5      32       9      39
 3J        7      12       7      13
 1C       13      42       8      35
 3B        3      33       6      37
 1M        6       8       4      25
 3C        6      27       5      23
 1O        8      26       4      23
 1J        6      18       6      16
 1T        1       3       3      19
 3E        2      12       4      14
 1G        1       7       4      11
 3D        7      22       4       9
 3G        2       5       1      10
 1L        2       4       2       8
 1I        1       3       1       7
 1P       16      65       0       3
 1F        0       1       0       3
 1S        1       2       1       2
 ;

 data csect95;
    set csection end=eof;
    if eof then call symput('ngroups', left( put( _n_, 4. )));
 run;

 proc means noprint data=csect95;
    var total95;
    output out=total sum=n;
 data total;
    set total end=eof;
    if eof then call symput('n', left( put( n, 4. )));
    put n;
 run;

 data csect95;
    set csection;
    keep id p95l p951 p95x p953 p95h p95n p95s p95m n95;
    p95x  = csect95 / total95;
    p95n  = total95;
    alpha = 0.01;
    g     = probit( 1 - ( alpha / ( 2 * &ngroups ) ) );
    p953  = p95x + g * sqrt( p95x * ( 1 - p95x ) / &n );
    p951  = p95x - g * sqrt( p95x * ( 1 - p95x ) / &n );

    z     = probit( 1 - ( alpha / 2 * &ngroups ));
    p95l  = p95x - z * sqrt( p95x * ( 1 - p95x ) / p95n );
    p953  = p95x + z * sqrt( p95x * ( 1 - p95x ) / p95n );

    /* assign dummy variables */
    p95h  = p953;
    p95l  = p951;
    p95s  = 0.001;
    p95m  = p95x;
 run;

 proc print;
 run;

 title 'C-Sections as a Proportion of Total'
       ' Deliveries in 1995';

 symbol v=none w=1;

 proc shewhart history=csect95;
    boxchart p95*id (n95) /
       nolimitslegend
       nolcl
       noucl
       xsymbol = 'Avg'
       turnhlabels
       blocklabtype = scaled
       nolegend
       npanel = &ngroups
       cframe = ligr
       cboxes = black
       haxis = axis1
       cboxfill = yellow;

 axis1 value = ( h=2.2 pct ) ;
 label p95x = 'Simultaneous 99% Intervals';

 run;