%newsas(sampdist);
title  'Generate sampling distribution of B0, B1';
*
*  X values are fixed, using the values from the Personality
*  test, Therapy example (X = personality test score).
*  For each X value, generate a set of sample Y's, according
*  to the model,
*           Y  =  14 + 1.2 * X + sigma * N(0,1)
*
*   where sigma = sqrt(MSE) = sqrt(80)
*         N(0,1)= a random normal deviate
*-----------------------------------------------------------;
options ls=80 ;
data XVALS;
     INPUT X @@;
     list;
CARDS;
26 24 22 33 27 36 30 38 30 34
DATA SAMPLES;
     set XVALS;                   *-- Read an X value;
     do SAMPLE=1 to 100;          *-- Generate Ys ;
         Y = 14 + 1.2 * X + SQRT(80)*NORMAL(16539);
         Y = round (Y);
         output;
     end;
run;
 
*---- Sort the data by SAMPLE;
proc sort data=SAMPLES;
     by SAMPLE;
Proc reg data = SAMPLES ;
     model Y = X;
     where sample=1;
     title2 'Analysis of the first sample';
 
/*------------------------------------------------------------*
 |       Fit a linear regression to each sample data,         |
 |  Save the estimates of B0, B1 in output data set, ESTIMATE |
 *------------------------------------------------------------*/
options nonotes;
Proc reg data = SAMPLES noprint
         outest= ESTIMATE;
     by SAMPLE;                   * Separate analysis for each sample;
     model Y = X;
     run;
options notes;
Title 'Regression estimates for each sample';
*---- Rename the slope and interecpt parameters;
data ESTIMATE;
     set ESTIMATE(rename=(intercep=b0 x=b1 _rmse_=sigma));
     drop  Y _TYPE_ _MODEL_ _DEPVAR_;
run;
data samples;
   merge samples estimate(keep=sample b0 b1);
   by sample;
   title 'First 4 samples and their estimates';
 
proc print data=SAMPLES;
     id SAMPLE b0 b1;
     by SAMPLE b0 notsorted b1 notsorted;
     where (SAMPLE<=4);
     title 'The first 4 random samples (Raw data)';
 
Proc print data=ESTIMATE;
     id sample;
     where (SAMPLE<=50);
     title 'Parameter estimates';
proc chart data=ESTIMATE;
     vbar B0 / midpoints = -45 to 55 by 10;
     vbar B1 / midpoints = -.4 to 2.8 by .4;
Proc plot;
     Plot B1 * B0 = SAMPLE /vref=1.2 href=14 vpos=40;
     title 'Joint sampling distribution of b0 & b1';
 
run;
 
*---- Determine empirical confidence limits;
proc means data=XVALS noprint;
     output out=XSUMRY  n=N css=CSS mean=XBAR;
 
data LIMITS;
     set ESTIMATE;
     retain N CSS XBAR;
     drop SIGMA N CSS XBAR XB2;
     drop T95 T975 beta0 beta1 sb0 sb1;
  *  drop c90b0 c90b1 b90b0 b90b1;
 
     beta0 = 14;             * True intercept;
     beta1 = 1.2;            *      & slope;
     T95 = 1.86 ;            * Critical t(.95,8) for alpha=.10 ;
     T975= 2.302;            * Bonferroni t for alpha=.10, g=2;
 
     if _N_=1 then SET XSUMRY;
     XB2  = XBAR * XBAR;
 
     MSE = SIGMA**2;         * Use sample MSE;
     SB0 = sqrt(MSE * ((1/N) + XB2/CSS));    * Std. error of B0;
     SB1 = sqrt(MSE / CSS);                  * Std. error of B1;
 
     c90b0 = t95 *sb0;     * Width of CI for B0;
     c90b1 = t95 *sb1;     * Width of CI for B1;
     B90b0 = t975*sb0;     * Width of Bonferroni interval for B0;
     B90b1 = t975*sb1;     * Width of Bonferroni interval for B1;
     low_b1= b1 - c90b1;
     high_b1=b1 + c90b1;
     low_b0= b0 - c90b0;
     high_b0=b0 + c90b0;
 
     label ERR_B0 = 'B0 outside CI ?'
           ERR_B1 = 'B1 outside CI ?'
           ERR_ANY= 'B1 or B0 outside CI ?'
           BERR_B0 = 'B0 outside Bonf. Int ?'
           BERR_B1 = 'B1 outside Bonf. Int ?'
           BERR_ANY= 'B1 or B0 outside Bonf. Int ?';
 
*--- Determine if sample limits include true value ;
     if low_b0 <= beta0 <= high_b0
              then ERR_B0 = 0;
              else ERR_B0 = 1;
     if low_b1 <= beta1 <= high_b1
              then ERR_B1 = 0;
              else ERR_B1 = 1;
     ERR_ANY = ERR_B0 OR ERR_B1;
 
*--- Determine if sample results are within Bonferroni Limits ;
     if (b0 - b90b0) <= beta0 <= (b0 + b90b0)
              then BERR_B0 = 0;
              else BERR_B0 = 1;
     if (b1 - b90b1) <= beta1 <= (b1 + b90b1)
              then BERR_B1 = 0;
              else BERR_B1 = 1;
     BERR_ANY = BERR_B0 OR BERR_B1;
     Title 'Empirical Confidence Limits';
proc freq;
     tables ERR_B0 ERR_B1 ERR_ANY /nocum;
     tables BERR_B0 BERR_B1 BERR_ANY /nocum;
 
proc timeplot data=limits;
   plot low_b1 ='<'
        high_b1='>'
        b1     ='*' / overlay hiloc npp ref=1.2
                     axis=-2 to 4 pos=0;
   id sample;
   run;