%newsas(normplt1);
/*-------------------------------------------------------------*
 |  Shows 4 ways to get a normal probability plot:             |
 |   (1) Direct calculation using PROBIT function in data step |
 |   (2) Using Proc Rank for Blom's approximation              |
 |   (3) Proc Univariate                                       |
 |   (4) normplot program  ** easiest **                       |
 *-------------------------------------------------------------*/
%include nwk(nwk01p20);   *-- Calculator maintenance data;
 
title 'Normal probability plot of residuals';
* Get residuals in output data set;
proc reg data=calc;
     model y = x;
     output out=results       /* name of output data set */
              r=resid         /* Residuals               */
              student=stdres  /* Std. residuals, e/sqrt(MSE) */
              pred=yhat;      /* Fitted values           */
run;
*---------------------------------------------------------------;
*(1) This data step and the Proc sort calculate the normal prob.;
*     values directly.  The PROC RANK step gets the same values;
*---------------------------------------------------------------;
proc sort data=results;
     by resid;
 
data quantile;
     set results;
     format p 8.4;
     p = (_N_ - .375) / 18.25;  * The proportion of residuals
                                    <= this one;
     Z  = Probit(p);            * Inverse normal function;
proc print;
     var x y yhat resid p z;
     title2 'Residuals and normal probability values (calculated)';
run;
*----------------------------------------------------------------*
| (2) Use proc RANK to get Normal quantiles. The BLOM measure is |
| defined to give the values used by Neter, Wasserman & Kutner:  |
|       E(resid[i]) = z ( (i-.375)/(n+.25) )                     |
| where z (p) is the inverse cumulative normal distribution      |
*----------------------------------------------------------------*;
proc rank normal=BLOM data=results out=quantile;
     var STDRES;          * The input variable;
     ranks Z;             * The output variable;
     Label Z = 'Expected Normal Quantile';
proc print data=quantile;
     var x y yhat resid stdres z;
     title2 'Residuals and normal probability values (Proc Rank)';
 
proc plot data=quantile;
     plot stdres*z = '*'          /* actual points    */
          z * z    = '+'          /* theoretical line */
     /  overlay vpos=30 hpos=60;
title2; run;
*----------------------------------------------------------------;
* (3) Proc Univariate also gives a normal probability plot if    ;
*  the PLOT option is used.  However, the size of the plot       ;
*  cannot be controlled (it depends on the pagesize and linesize);
 
proc univariate PLOT;
     var stdres;
run;
*--------------------------------------------------------------;
* (4) The NORMPLOT program provides a normal probability plot  ;
*  with the theoretical line, and a detrended version, rotated ;
*  so the theoretical line is flat.                            ;
*--------------------------------------------------------------;
*include macros(normplot);
%normplot(data=results,var=resid);