/* arsenic2.sas */
options ls=72 ps=60;
data arsenic;
   infile 'arsenic.dat';
   input conc dist;
   logconc = log(conc);
proc plot data=arsenic;
   plot logconc*dist;
title 'plot of logconc vs dist';   
proc glm data=arsenic;
   model logconc = dist / ss1 ss2 p clm clparm ;
      output out=new predicted=yhat residual=resid lcl=clpl
             ucl=clpu lclm=clml uclm=clmu;
title 'regression of logconc on dist';
proc print data=new;
   var conc dist clml clmu clpl clpu yhat resid;
proc plot data=new;
   plot resid*dist / vref=0;
proc univariate data=new plot normal;
   var resid;
title 'descriptive statistics for residuals';
proc rank data=new normal=blom out=new1;
   var resid;
   ranks nscore;
proc plot data=new1;
   plot nscore*resid;
title 'normal probability plot for residuals';
/* get estimate and prediction for dist = 20 */
data dist20;
   input conc dist;
   logconc = log(conc);   
cards;
. 20
;
data combine;
   set arsenic dist20;   
proc glm data=combine;
   model logconc = dist / ss1 ss2 p clm clparm ;
      output out=new1a predicted=yhat residual=resid lcl=clpl
             ucl=clpu lclm=clml uclm=clmu;
title 'regression of logconc on dist (2)';
proc print data=new1a;
   var conc dist clml clmu clpl clpu yhat resid;
/* backtransform */
data back;
   set new1a;
   clmlback = exp(clml);
   clmuback = exp(clmu);
   clplback = exp(clpl);
   clpuback = exp(clpu);
   yhatback = exp(yhat);
proc print data=back;
   var conc dist clmlback clmuback clplback clpuback yhatback; 
title 'backtransformed values';
/* add data and form plots */
data extra;
   input conc dist;
   logconc = log(conc);
cards;
. 2
. 4
. 6
. 8
. 10
. 12
. 14
. 16
. 18
. 20
. 22
. 24
. 26
. 28
. 30
. 32
. 34
. 36
;
data combine2;
   set arsenic extra; 
proc glm data=combine2 noprint;
   model logconc = dist;
      output out=new2a predicted=yhat residual=resid;
title 'plots for regression of logconc on dist';
data moreplot;
   set new2a;
   expyhat = exp(yhat);
   code = '*';
   if conc = . then code = '+';
   yhat2 = logconc;
   if conc = . then yhat2 = yhat;
   expyhat2 = conc;
   if conc = . then expyhat2 = expyhat;
proc print data=moreplot;
   var conc dist yhat expyhat yhat2 expyhat2 code;
proc plot data=moreplot;
   plot yhat2*dist=code;
   plot expyhat2*dist=code;
run;